{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-VRZnOBNA6o7"
   },
   "source": [
    "# Regression Discontinuity\n",
    "This notebook illustrates the use of Regression Discontinuity in an empirical study. We analyze the effect of the antipoverty program *Progresa/Oportunidades* on the consumption behavior of families in Mexico in the early 2000s.\n",
    "\n",
    "The program was intended for families in extreme poverty and included financial incentives for participation in measures that improved the family's health, nutrition and children's education. The effect of this program is a widely studied problem in social and economic sciences and, according to the WHO, was a very successful measure in terms of reducing extreme poverty in Mexico.\n",
    "\n",
    "Eligibility for the program was determined based on a pre-intervention household poverty-index. Individuals above a certain threshold received the treatment (participation in the program) while individuals below the threshold were excluded and recorded as a control group. All observations above the threshold participated in the program, which makes the analysis fall into the standard (sharp) regression discontinuity design.\n",
    "\n",
    "First, we need to install and load some packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "1Yr5aL2yAgYN"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defaulting to user installation because normal site-packages is not writeable\n",
      "Requirement already satisfied: rdd in c:\\users\\vasilis\\appdata\\roaming\\python\\python311\\site-packages (0.0.3)\n",
      "Requirement already satisfied: rdrobust in c:\\users\\vasilis\\appdata\\roaming\\python\\python311\\site-packages (1.2.1)\n",
      "Requirement already satisfied: pandas in c:\\users\\vasilis\\appdata\\roaming\\python\\python311\\site-packages (from rdd) (2.2.2)\n",
      "Requirement already satisfied: numpy in c:\\programdata\\anaconda3\\lib\\site-packages (from rdd) (1.24.4)\n",
      "Requirement already satisfied: statsmodels in c:\\programdata\\anaconda3\\lib\\site-packages (from rdd) (0.14.0)\n",
      "Requirement already satisfied: scipy in c:\\programdata\\anaconda3\\lib\\site-packages (from rdrobust) (1.11.1)\n",
      "Requirement already satisfied: scikit-learn>=1.2.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from rdrobust) (1.2.2)\n",
      "Requirement already satisfied: plotnine in c:\\users\\vasilis\\appdata\\roaming\\python\\python311\\site-packages (from rdrobust) (0.13.6)\n",
      "Requirement already satisfied: python-dateutil>=2.8.2 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas->rdd) (2.8.2)\n",
      "Requirement already satisfied: pytz>=2020.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas->rdd) (2023.3.post1)\n",
      "Requirement already satisfied: tzdata>=2022.7 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas->rdd) (2023.3)\n",
      "Requirement already satisfied: joblib>=1.1.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from scikit-learn>=1.2.0->rdrobust) (1.2.0)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from scikit-learn>=1.2.0->rdrobust) (2.2.0)\n",
      "Requirement already satisfied: matplotlib>=3.7.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from plotnine->rdrobust) (3.7.2)\n",
      "Requirement already satisfied: mizani~=0.11.0 in c:\\users\\vasilis\\appdata\\roaming\\python\\python311\\site-packages (from plotnine->rdrobust) (0.11.4)\n",
      "Requirement already satisfied: patsy>=0.5.2 in c:\\programdata\\anaconda3\\lib\\site-packages (from statsmodels->rdd) (0.5.3)\n",
      "Requirement already satisfied: packaging>=21.3 in c:\\programdata\\anaconda3\\lib\\site-packages (from statsmodels->rdd) (23.1)\n",
      "Requirement already satisfied: contourpy>=1.0.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from matplotlib>=3.7.0->plotnine->rdrobust) (1.0.5)\n",
      "Requirement already satisfied: cycler>=0.10 in c:\\programdata\\anaconda3\\lib\\site-packages (from matplotlib>=3.7.0->plotnine->rdrobust) (0.11.0)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from matplotlib>=3.7.0->plotnine->rdrobust) (4.25.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from matplotlib>=3.7.0->plotnine->rdrobust) (1.4.4)\n",
      "Requirement already satisfied: pillow>=6.2.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from matplotlib>=3.7.0->plotnine->rdrobust) (9.4.0)\n",
      "Requirement already satisfied: pyparsing<3.1,>=2.3.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from matplotlib>=3.7.0->plotnine->rdrobust) (3.0.9)\n",
      "Requirement already satisfied: six in c:\\programdata\\anaconda3\\lib\\site-packages (from patsy>=0.5.2->statsmodels->rdd) (1.16.0)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\vasilis\\AppData\\Roaming\\Python\\Python311\\site-packages\\pandas\\core\\arrays\\masked.py:60: UserWarning: Pandas requires version '1.3.6' or newer of 'bottleneck' (version '1.3.5' currently installed).\n",
      "  from pandas.core import (\n"
     ]
    }
   ],
   "source": [
    "!pip install rdd rdrobust\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LinearRegression, LassoCV\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from lightgbm import LGBMRegressor\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import patsy\n",
    "import math\n",
    "from rdd.rdd import optimal_bandwidth\n",
    "from rdrobust import rdrobust"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GH0wFmHSxnen"
   },
   "source": [
    "We use a dataset assembled by [Calonico et al. (2014)](https://rdpackages.github.io/references/Calonico-Cattaneo-Titiunik_2014_ECMA--Supplemental.pdf) and follow the analysis in [Noack et al. (2023)](https://arxiv.org/pdf/2107.07942.pdf).\n",
    "\n",
    "First, we open the data and remove any observations that have NaN values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "Rzbv0XXCxxJt"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of Data:\n",
      "(1944, 27)\n",
      "Variable Names:\n",
      "Index(['hhpiso', 'hhrooms', 'hhwater', 'hhwaterin', 'hhbano', 'hhownhouse',\n",
      "       'hhsize', 'hhelect', 'clus', 'headmale', 'headage', 'heademp',\n",
      "       'wifeage', 'wifeeduc', 'headeduc', 'child_0to5', 'boy_0to5',\n",
      "       'pov_index', 'conspcfood_t0', 'conspcfood_t1', 'conspcfood_t2',\n",
      "       'conspcnonfood_t0', 'conspcnonfood_t1', 'conspcnonfood_t2', 'conspc_t0',\n",
      "       'conspc_t1', 'conspc_t2'],\n",
      "      dtype='object')\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>hhpiso</th>\n",
       "      <th>hhrooms</th>\n",
       "      <th>hhwater</th>\n",
       "      <th>hhwaterin</th>\n",
       "      <th>hhbano</th>\n",
       "      <th>hhownhouse</th>\n",
       "      <th>hhsize</th>\n",
       "      <th>hhelect</th>\n",
       "      <th>clus</th>\n",
       "      <th>headmale</th>\n",
       "      <th>...</th>\n",
       "      <th>pov_index</th>\n",
       "      <th>conspcfood_t0</th>\n",
       "      <th>conspcfood_t1</th>\n",
       "      <th>conspcfood_t2</th>\n",
       "      <th>conspcnonfood_t0</th>\n",
       "      <th>conspcnonfood_t1</th>\n",
       "      <th>conspcnonfood_t2</th>\n",
       "      <th>conspc_t0</th>\n",
       "      <th>conspc_t1</th>\n",
       "      <th>conspc_t2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>1.5410</td>\n",
       "      <td>260.223999</td>\n",
       "      <td>374.071991</td>\n",
       "      <td>379.208008</td>\n",
       "      <td>106.186668</td>\n",
       "      <td>183.149323</td>\n",
       "      <td>57.127998</td>\n",
       "      <td>366.410675</td>\n",
       "      <td>557.221313</td>\n",
       "      <td>436.335999</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.6615</td>\n",
       "      <td>265.359985</td>\n",
       "      <td>574.375977</td>\n",
       "      <td>195.167999</td>\n",
       "      <td>87.266663</td>\n",
       "      <td>160.013336</td>\n",
       "      <td>112.136002</td>\n",
       "      <td>352.626648</td>\n",
       "      <td>734.389282</td>\n",
       "      <td>307.304016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.3240</td>\n",
       "      <td>604.549988</td>\n",
       "      <td>433.350006</td>\n",
       "      <td>367.010010</td>\n",
       "      <td>348.813324</td>\n",
       "      <td>263.083344</td>\n",
       "      <td>123.323334</td>\n",
       "      <td>953.363281</td>\n",
       "      <td>696.433350</td>\n",
       "      <td>490.333344</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>1.1500</td>\n",
       "      <td>233.259995</td>\n",
       "      <td>281.766663</td>\n",
       "      <td>260.468567</td>\n",
       "      <td>100.877777</td>\n",
       "      <td>140.233337</td>\n",
       "      <td>36.325714</td>\n",
       "      <td>334.137756</td>\n",
       "      <td>422.000000</td>\n",
       "      <td>296.794281</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.8260</td>\n",
       "      <td>280.339996</td>\n",
       "      <td>185.110001</td>\n",
       "      <td>364.869995</td>\n",
       "      <td>158.416656</td>\n",
       "      <td>201.183334</td>\n",
       "      <td>68.000000</td>\n",
       "      <td>438.756653</td>\n",
       "      <td>386.293335</td>\n",
       "      <td>432.869995</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 27 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   hhpiso  hhrooms  hhwater  hhwaterin  hhbano  hhownhouse  hhsize  hhelect  \\\n",
       "1     1.0      1.0      1.0        1.0     1.0         0.0       5        1   \n",
       "2     1.0      1.0      1.0        0.0     0.0         1.0       5        1   \n",
       "3     1.0      1.0      0.0        0.0     1.0         1.0       4        1   \n",
       "4     0.0      1.0      0.0        0.0     1.0         0.0       6        1   \n",
       "6     1.0      1.0      1.0        0.0     1.0         1.0       4        1   \n",
       "\n",
       "   clus  headmale  ...  pov_index  conspcfood_t0  conspcfood_t1  \\\n",
       "1     1         1  ...     1.5410     260.223999     374.071991   \n",
       "2     1         1  ...    -0.6615     265.359985     574.375977   \n",
       "3     1         1  ...    -0.3240     604.549988     433.350006   \n",
       "4     1         1  ...     1.1500     233.259995     281.766663   \n",
       "6     1         1  ...    -0.8260     280.339996     185.110001   \n",
       "\n",
       "   conspcfood_t2  conspcnonfood_t0  conspcnonfood_t1  conspcnonfood_t2  \\\n",
       "1     379.208008        106.186668        183.149323         57.127998   \n",
       "2     195.167999         87.266663        160.013336        112.136002   \n",
       "3     367.010010        348.813324        263.083344        123.323334   \n",
       "4     260.468567        100.877777        140.233337         36.325714   \n",
       "6     364.869995        158.416656        201.183334         68.000000   \n",
       "\n",
       "    conspc_t0   conspc_t1   conspc_t2  \n",
       "1  366.410675  557.221313  436.335999  \n",
       "2  352.626648  734.389282  307.304016  \n",
       "3  953.363281  696.433350  490.333344  \n",
       "4  334.137756  422.000000  296.794281  \n",
       "6  438.756653  386.293335  432.869995  \n",
       "\n",
       "[5 rows x 27 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"https://raw.githubusercontent.com/CausalAIBook/MetricsMLNotebooks/main/data/progresa.csv\",\n",
    "                 index_col=0)\n",
    "df = df.dropna()\n",
    "df.rename(columns={\"index\": \"pov_index\"}, inplace=True)\n",
    "print(\"Shape of Data:\")\n",
    "print(df.shape)\n",
    "print(\"Variable Names:\")\n",
    "print(df.columns)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vGbvqQmpmoqV"
   },
   "source": [
    "The data set contains 1,944 observations for which full covariate information of 27 variables is available.\n",
    "\n",
    "We want to measure the local average treatment effect of program participation on four outcome variables. The outcome variables are food and non-food consumption of the recorded families at two points in time, one year and two years after the implementation of the program.\n",
    "\n",
    "The baseline covariates, recorded prior to program implementation, include the household's size; household head's age, sex, years of education and employment status; spouse's age and years of education; number of children not older than five years and their sex, and physical characteristics of the house: whether the house has cement floors, water connection, water connection inside the house, a bathroom, electricity, number of rooms, pre-intervention consumption, and an identifier of the urban locality in which the house is located.\n",
    "\n",
    "The data fits to the pattern of a sharp RD design, namely, all individuals that were below the cut-off index received no intervention, and all individuals above the cut-off were eligible to join the *progresa* program and thus participated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9yvX75wy98g9"
   },
   "source": [
    "## Estimation without Covariates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bCueRzpuqNXn"
   },
   "source": [
    "First, we will perform a very simple RD estimation with a weighted linear regression. We use a triangular kernel, which assigns weights to observations based on their distance from the cutoff point. The weights decrease linearly as the distance from the cutoff point increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "1LAMZP540pLn"
   },
   "outputs": [],
   "source": [
    "def triangular_kernel(index, h):\n",
    "    weights = 1 - np.abs(index) / h\n",
    "    weights[weights < 0] = 0\n",
    "    return weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "N-I-EBps0ubO"
   },
   "source": [
    "The parameter `h` is the bandwidth that controls the range of observations that receive non-zero weights. We use the `IKbandwidth` function from the `rdd` package that implements the *Imbens-Kalyanaraman* method. Another standard approach would be to use the standard deviation of `index`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "bFuzAouP04lO"
   },
   "outputs": [],
   "source": [
    "h = optimal_bandwidth(X=df.pov_index, Y=df.conspcfood_t1, cut=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "J9kU7tQ207A3"
   },
   "source": [
    "We use the triangular kernel function to calculate weights for each observation. After that, we can fit two seperate linear regressions for both treatment and control groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "cjc7f7F6qM36"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>LinearRegression()</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">LinearRegression</label><div class=\"sk-toggleable__content\"><pre>LinearRegression()</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights = triangular_kernel(df.pov_index, h)\n",
    "model_treated, model_control = LinearRegression(), LinearRegression()\n",
    "model_treated.fit(y=df.loc[df.pov_index > 0, \"conspcfood_t1\"].values.reshape(-1, 1),\n",
    "                  X=df.loc[df.pov_index > 0, \"pov_index\"].values.reshape(-1, 1),\n",
    "                  sample_weight=weights[df.pov_index > 0])\n",
    "model_control.fit(y=df.loc[df.pov_index < 0, \"conspcfood_t1\"].values.reshape(-1, 1),\n",
    "                  X=df.loc[df.pov_index < 0, \"pov_index\"].values.reshape(-1, 1),\n",
    "                  sample_weight=weights[df.pov_index < 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MC5vPB-I1jeH"
   },
   "source": [
    "The treatment effect at the cutoff point is estimated as the difference between the predictions of the two models at the cutoff point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "279my1C8o9a3"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-22.105321594809766"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cutoff = 0\n",
    "treatment_effect = model_treated.predict(np.array([cutoff]).reshape(-1, 1))\n",
    "treatment_effect -= model_control.predict(np.array([cutoff]).reshape(-1, 1))\n",
    "treatment_effect[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uW6PYdz-BESB"
   },
   "source": [
    "We estimate that the participation in the program reduced food consumption by $22.1$ units in the year following the intervention.The following plot visualizes the two weighted regressions at the cut-off for the last outcome variable (for food consumption in `t1`). We can clearly see the \"jump\" at the cut-off, which is our LATE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "uvONPb44A56z"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAGdCAYAAAAWp6lMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABXAUlEQVR4nO3deXhTVd4H8G+6JN1DF2gplBYGSoWyIw2gLLKJIiAywMhL0WFQlEVkaUBHARVpq4g4I4zDMCAdtS4I47hUFikutWFHoEAFyyatKEJLS5tu5/0jNDQkLU2b5CY338/z3Iebc2/uPadJyC9nVQghBIiIiIhkzEPqDBARERHZGwMeIiIikj0GPERERCR7DHiIiIhI9hjwEBERkewx4CEiIiLZY8BDREREsseAh4iIiGTPS+oMNEZ1dTUuXryIwMBAKBQKqbNDREREDSCEwLVr1xAZGQkPD8fWubhkwHPx4kVERUVJnQ0iIiJqhPPnz6N169YOvadLBjyBgYEADH+woKAgiXNDRDZTUgJERhr2L14E/P2lzQ8R2VRRURGioqKM3+OO5JIBT00zVlBQEAMeIjnx9Ly5HxTEgIdIpqTojsJOy0RERCR7LlnDQ0Qy5eUFTJ16c5+IyEb4PwoROQ+VCti4UepcEJEMMeAhIpKYEAKVlZWoqqqSOitkB4cPH8aZM2cQExODbt26SZ0du/P29oZn7f54ToIBDxE5DyGA69cN+35+gBvMs1VeXo78/Hxcryk3ycqVK1dQVFQEtVqNK1eu4MCBAwgODpY6W3alUCjQunVrBAQESJ0VEwx4iMh5XL8O1PwnWVws+1Fa1dXVyMvLg6enJyIjI6FUKjmZqoyUlJSgpKQEYWFhJuktWrSAv0zf20II/Prrr7hw4QI6dOjgVDU9DHiIiCRSXl6O6upqREVFwc/PT+rskI2VlJTUeczHx8eBOXGs5s2b48yZM6ioqHCqgIfD0omIJOboKfbJMVQqlVXpcuGstZT8lBEREdlBQEAAIiIiTNIiIiKcrm+Lu2DAQ0REbiEzMxMKhQJXr1512D1bt26NuLg4tG3bFnFxcQ5fP6qGFGV3Ngx4iMjt6HQ6pKWlQafTSZ0Vl1VQUIDZs2ejXbt2UKlUiIqKwgMPPICdO3fa9D6DBg3C3LlzbXrN+sTExEChUEChUMDX1xdxcXF45ZVXIIRo9DUDAgIQGhoqac1Ov379kJ+fD7VaLVkepMZOy0TkVrRaLVJTU42Pk5KSkJKSImGOXM+ZM2fQv39/NGvWDKmpqejatSsqKirw5ZdfYubMmThx4oRD8yOEQFVVFbxsNDv3Cy+8gOnTp6OsrAw7duzAE088gaCgIDz++OM2ub4lFRUV8Pb2ttv1lUqlWfOau2END7kV/rJ3cp6ewPjxhs0Oozt0Op1JsAMAqampTvN+qK4WuFysl2yrrm5YLcaTTz4JhUKBPXv2YPz48YiNjUXnzp0xb948ZGdnG887d+4cxowZg4CAAAQFBWHChAn45ZdfjMeXLl2K7t27Iy0tDTExMVCr1Zg0aRKuXbsGAHjkkUewe/durF692ljrcubMGWPzzJdffonevXtDpVLhm2++gV6vx5w5c9CiRQv4+Pjgrrvuwt69e61+HQIDAxEREYGYmBj85S9/QdeuXbFt2zbj8fLyciQlJaFVq1bw9/dHQkICMjMzTa6xbt064+i7Bx98EK+99hqaNWtmVvZ///vfxloyIQQKCwvx2GOPoUWLFggKCsI999yDw4cPG593+PBhDB48GIGBgQgKCkKvXr2wb98+AMDZs2fxwAMPIDg4GP7+/ujcuTM+//xzAJabtDZv3ozOnTtDpVIhJiYGK1euNClDTEwMXn75Zfz5z39GYGAg2rRpg3/+859W/z2dBWt4yG3wl70L8PEBPvzQbpfPzc2tMz0hIcFu922oK9fL0eulHZLdf/9fhyI0oP4RRL///jsyMjKwfPlyi3PJ1HypCyEwduxY+Pv7Y/fu3aisrMSTTz6JiRMnmgQHp0+fxtatW/Hpp5/iypUrmDBhApKTk7F8+XKsXr0aubm5iI+PxwsvvADg5pBnwPAZfvXVV9GuXTs0a9YMSUlJ2Lx5M95++21ER0cjNTUVI0aMwKlTpxASEmL130MIgd27d+P48ePo0KGDMf3RRx/FmTNnkJ6ejsjISGzZsgX33nsvjhw5gg4dOuC7777DjBkzkJKSgtGjR2PHjh147rnnzK5/6tQpfPDBB9i8ebNx+Pb999+PkJAQfP7551Cr1XjrrbcwZMgQ5ObmIiQkBJMnT0aPHj2wdu1aeHp64tChQ8aaoZkzZ6K8vBxff/01/P39kZOTU2cz2v79+zFhwgQsXboUEydORFZWFp588kmEhobikUceMZ63cuVKvPjii3jmmWfw0Ucf4YknnsCAAQMQFxdn9d9Tagx4yC3U9ct+3Lhxkn7R6XQ65ObmIjY21iwf9R2jxomNjbUqncydOnUKQojbfuHt2LEDP/zwA/Ly8hAVFQUASEtLQ+fOnbF3717ceeedAAyTL27cuBGBgYEAgClTpmDnzp1Yvnw51Go1lEol/Pz8LDbHvPDCCxg2bBgAw5w3a9euxcaNGzFy5EgAhlqW7du3Y/369Vi4cGGDy6jVavHXv/4V5eXlqKiogI+PD+bMmQPAEKC99957uHDhAiIjIwEACxYsQEZGBjZs2ICXX34Zf/vb3zBy5EgsWLAAgOH9lZWVhU8//dTkPuXl5UhLS0Pz5s0BAF999RWOHDmCS5cuGYeuv/rqq9i6dSs++ugjPPbYYzh37hwWLlxo/PvXDsTOnTuHhx56CF26dAEAtGvXrs4yvvbaaxgyZIgxEIuNjUVOTg5eeeUVk4Dnvvvuw5NPPmn8u6xatQqZmZkuGfCwSYvcQn2/7KWi1Wqh0WiQmJgIjUYDrVbboGPUeAkJCUhKSjJJ02q1DCitUNN593ZzrRw/fhxRUVHGYAcAOnXqhGbNmuH48ePGtJiYGGOwAwAtW7bEpUuXGpSX3r17G/dPnz6NiooK9O/f35jm7e2NPn36mNyvIRYuXIhDhw5h9+7dGDx4MJ599ln069cPAHDgwAEIIRAbG4uAgADjtnv3bpw+fRoAcPLkSfTp08fkmrc+BoDo6GhjsAMYal2Ki4uNHZxrtry8POO1582bh7/85S8YOnQokpOTjekAMGfOHLz00kvo378/lixZgh9++KHOMh4/ftzkbwUA/fv3x48//miyplvXrl2N+wqFAhEREQ1+fZwNa3jILTjbL/v6apxq9i0dk/0Xc0mJ3ZeWSElJwbhx41h71kgdOnSAQqHA8ePHMXbs2DrPE0JYDIpuTb+1o65CoUB1dXWD8lK7Sa2uQKyufNQnLCwM7du3R/v27bF582a0b98eGo0GQ4cORXV1NTw9PbF//36zWYRrmo8s3dPSKK9bmwSrq6vRsmVLs/5AwM2mwqVLl+Lhhx/GZ599hi+++AJLlixBeno6HnzwQfzlL3/BiBEj8Nlnn2Hbtm1YsWIFVq5cidmzZ5tdr6F5bMrr42wY8JBbqPllXzuQkPKXfWNqnJyln4kcJCQkOOXfMthPif1/HSrp/W8nJCQEI0aMwJtvvok5c+aYfWlfvXoVzZo1Q6dOnXDu3DmcP3/eWMuTk5ODwsJC3HHHHQ3Ok1KpbNAq8u3bt4dSqcS3336Lhx9+GIBh5NO+ffuaNKw9ODgYs2fPxoIFC3Dw4EH06NEDVVVVuHTpEu6++26Lz4mLi8OePXtM0mo6FtenZ8+eKCgogJeXF2JiYuo8LzY2FrGxsXj66afxpz/9CRs2bMCDDz4IAIiKisKMGTMwY8YMLF68GOvWrbMY8HTq1AnffvutSVpWVhZiY2OdajkIW2LAQ27DmX7ZN6bGif1M5M/DQ3HbTsPOYM2aNejXrx/69OmDF154AV27dkVlZSW2b9+OtWvX4vjx4xg6dCi6du2KyZMn4/XXXzd2Wh44cKBJU9TtxMTEQKfT4cyZMwgICKiz87G/vz+eeOIJLFy4ECEhIWjTpg1SU1Nx/fp1TJs2rUnlnTlzJlJSUrB582aMHz8ekydPRmJiIlauXIkePXrgt99+w1dffYUuXbrgvvvuw+zZszFgwAC89tpreOCBB/DVV1/hiy++uG1N09ChQ9G3b1+MHTsWKSkp6NixIy5evIjPP/8cY8eORefOnbFw4UKMHz8ebdu2xYULF7B371489NBDAIC5c+di5MiRiI2NxZUrV/DVV1/VGVzOnz8fd955J1588UVMnDgR33//Pf7+979jzZo1TfpbOTXhggoLCwUAUVhYKHVWiBotKSlJADBuWq22QcdkrbhYCMCwFRdLnRu7Ky0tFTk5OaK0tFTqrFjt4sWLYubMmSI6OloolUrRqlUrMXr0aLFr1y7jOWfPnhWjR48W/v7+IjAwUPzxj38UBQUFxuNLliwR3bp1M7nuqlWrRHR0tPHxyZMnhUajEb6+vgKAyMvLE7t27RIAxJUrV0yeW1paKmbPni3CwsKESqUS/fv3F3v27DEer+t5tUVHR4tVq1aZpU+fPl107txZVFVVifLycvH888+LmJgY4e3tLSIiIsSDDz4ofvjhB+P5//znP0WrVq2Er6+vGDt2rHjppZdEREREvWUXQoiioiIxe/ZsERkZKby9vUVUVJSYPHmyOHfunNDr9WLSpEkiKipKKJVKERkZKWbNmmV8/8yaNUv84Q9/ECqVSjRv3lxMmTJF/Pbbb3WW/aOPPhKdOnUS3t7eok2bNuKVV1657d+iW7duYsmSJXX+/YSo/30t5fe3QogmTB8pkaKiIqjVahQWFiIoKEjq7BA1Gkdp3cIBfXicSVlZGfLy8tC2bVtZr55NwPTp03HixAl88803UmfF7up7X0v5/c0mLSIJ1deXxFn7mRDR7b366qsYNmwY/P398cUXX+Dtt9+Wd3ORC2DAQ0REZGN79uxBamoqrl27hnbt2uGNN97AX/7yF6mz5dYY8BCR8/D0BO677+Y+kYv64IMPpM4C3YIBDxE5Dx8f4LPPpM4FEckQZ1omIiIi2WPAQ0RERLLHgIeInEdJiWEour+/YZ+IyEbYh4eInMv161LngIhkyKoanrVr16Jr164ICgpCUFAQ+vbtiy+++MLiuY8//jgUCgVef/11k3S9Xo/Zs2cjLCwM/v7+GD16NC5cuNDoAhARERHdjlUBT+vWrZGcnIx9+/Zh3759uOeeezBmzBgcO3bM5LytW7dCp9MhMjLS7Bpz587Fli1bkJ6ejm+//RbFxcUYNWpUgxaHIyJ50Ol0SEtLg06nkzor5IaWLl2K7t27S50NSdi67MXFxbh8+TKKi4ttdk17sSrgeeCBB3DfffcZV2pdvnw5AgICkJ2dbTzn559/xqxZs/DOO++YLStfWFiI9evXY+XKlRg6dCh69OiB//znPzhy5Ah27NhhmxIRkVPTarXQaDRITEyERqOBVquVOktkBYVCUe/2yCOP2OW+jgxSzpw5Y1ImtVoNjUaD//3vfw65vz0tWLAAO3futMm1Lly4gBMnTiAvLw8nTpxw+taaRndarqqqQnp6OkpKStC3b18AQHV1NaZMmYKFCxeic+fOZs/Zv38/KioqMHz4cGNaZGQk4uPjkZWV1disEJGL0Ol0SE1NNUlLTU1lTY8Lyc/PN26vv/46goKCTNJWr15tcn5FRYVEOW26HTt2ID8/HzqdDn369MFDDz2Eo0eP2vWe5eXldr1+QEAAQkNDm3yd4uJiFBQUmKQVFBQ4dU2P1QHPkSNHEBAQAJVKhRkzZmDLli3o1KkTACAlJQVeXl6YM2eOxecWFBRAqVQiODjYJD08PNzsD1ebXq9HUVGRyUZEric3N9eqdLdTXQ2U/CbdVl192yxGREQYN7VaDYVCYXxcVlaGZs2a4YMPPsCgQYPg4+OD//znPwCADRs24I477oCPjw/i4uLM1pXSarWIjY2Fn58f2rVrh+eee84YLG3cuBHLli3D4cOHjbUuGzduBGBoOXjsscfQokULBAUF4Z577sHhw4dNrp2cnIzw8HAEBgZi2rRpKCsra9DLERoaioiICMTFxWH58uWoqKjArl27jMd//vlnTJw4EcHBwQgNDcWYMWNw5swZ4/HKykrMmTMHzZo1Q2hoKLRaLaZOnYqxY8cazxk0aBBmzZqFefPmISwsDMOGDQMA5OTk4L777kNAQADCw8MxZcoU/Pbbb8bnffTRR+jSpQt8fX0RGhqKoUOHouTGyMbMzEz06dMH/v7+aNasGfr374+zZ88CMK8pq66uxgsvvIDWrVtDpVKhe/fuyMjIMB6vqe36+OOPMXjwYPj5+aFbt2749ttvLf7N9Hp9g/62UrB6lFbHjh1x6NAhXL16FZs3b8bUqVOxe/dulJaWYvXq1Thw4AAUCoVV1xRC1PucFStWYNmyZdZmlQiAm6467qRiY2PrTNfpdDh99ChG9ehhWEXZw3VmzbDZe6z0d+CVP9guY9ZaeBrwD2vyZbRaLVauXIkNGzZApVJh3bp1WLJkCf7+97+jR48eOHjwIKZPnw5/f39MnToVABAYGIiNGzciMjISR44cwfTp0xEYGIikpCRMnDgRR48eRUZGhrH7g1qthhAC999/P0JCQvD5559DrVbjrbfewpAhQ5Cbm4uQkBB88MEHWLJkCd58803cfffdSEtLwxtvvIF27do1uDwVFRVYt24dABi7aly/fh2DBw/G3Xffja+//hpeXl546aWXcO+99+KHH36AUqlESkoK3nnnHWOwt3r1amzduhWDBw82uf7bb7+NJ554At999x2EEMjPz8fAgQMxffp0vPbaaygtLYVWq8WECRPw1VdfIT8/H3/605+QmpqKBx98ENeuXcM333wDIQQqKysxduxYTJ8+He+99x7Ky8uxZ8+eOr9jV69ejZUrV+Ktt95Cjx498O9//xujR4/GsWPH0KFDB+N5zz77LF599VV06NABzz77LB5//HG8//778PIyDSNUKlWD/64OJ5poyJAh4rHHHhOrVq0SCoVCeHp6GjcAwsPDQ0RHRwshhNi5c6cAIH7//XeTa3Tt2lU8//zzdd6jrKxMFBYWGrfz588LAKKwsLCp2SeZS0pKEgCMW1JSktRZcnu3viZardalX6em5L20tFTk5OSI0tJSQ0Lxr0IsCZJuK/7VqrJv2LBBqNVq4+O8vDwBQLz++usm50VFRYl3333XJO3FF18Uffv2rfPaqampolevXsbHS5YsEd26dTM5Z+fOnSIoKEiUlZWZpP/hD38Qb731lhBCiL59+4oZM2aYHE9ISDC7Vm015fD19RX+/v7Cw8NDABAxMTHi8uXLQggh1q9fLzp27Ciqq6uNz9Pr9cLX11d8+eWXQgghwsPDxSuvvGI8XllZKdq0aSPGjBljTBs4cKDo3r27yf2fe+45MXz4cJO0mu+9kydPiv379wsA4syZM2Z5v3z5sgAgMjMzLZbt1r9jZGSkWL58uck5d955p3jyySdN/hb/+te/jMePHTsmAIhdu3aJvXv3Grfz588LISy8r2spLCyU7Pu7yfPwCCGg1+sxZcoUDB061OTYiBEjMGXKFDz66KMAgF69esHb2xvbt2/HhAkTABjag48ePWrWrl+bSqVy7qiRnFJd/UXGjRvHmh4JpaSkYNy4ccYaEQDQaDQm57jK68T3mGW9e/c27v/66684f/48pk2bhunTpxvTKysroVarjY8/+ugjvP766zh16hSKi4tRWVlpqOmrx/79+1FcXGzWJ6W0tBSnT58GABw/fhwzZswwOd63b1+Tpqm6vP/++4iLi0Nubi7mzp2Lf/zjHwgJCTHe+9SpUwgMDDR5TllZGU6fPo3CwkL88ssv6NOnj/GYp6cnevXqhepbmg5r/71qrr1r1y4EBASY5en06dMYPnw4hgwZgi5dumDEiBEYPnw4xo8fj+DgYISEhOCRRx7BiBEjMGzYMAwdOhQTJkxAy5Ytza5VVFSEixcvon///ibp/fv3N2sW7Nq1q3G/5loeHh6Ii4uDXq+HSqWymF9nYlXA88wzz2DkyJGIiorCtWvXkJ6ejszMTGRkZCA0NNTsTeft7Y2IiAh07NgRgKEKctq0aZg/fz5CQ0MREhKCBQsWoEuXLmbBElFT1ddfxJ2/jJxBQkKC8TVIS0uzeI4rvE58j1nm7+9v3K/5cl+3bp3Z38TT0xMAkJ2djUmTJmHZsmUYMWIE1Go10tPTsXLlynrvU11djZYtWyIzM9PsWLNmzZpWCABRUVHo0KEDOnTogICAADz00EPIyclBixYtUF1djV69euGdd94xe17z5s2N+7c2JQkhzM6v/fcCDOV64IEHkJKSYnZuy5Yt4enpie3btyMrKwvbtm3D3/72Nzz77LPQ6XRo27YtNmzYgDlz5iAjIwPvv/8+/vrXv2L79u1mPyzqy+OtabVHXdccq66uRkBAgNMHOjWsCnh++eUXTJkyBfn5+VCr1ejatSsyMjKMnawaYtWqVfDy8sKECRNQWlqKIUOGYOPGjcY3PpGt1NdfhJxH7dfDD8CZG/tnoqKkyI5VbP4e8w0x9KORim+IzS8ZHh6OVq1a4aeffsLkyZMtnvPdd98hOjoazz77rDGtppNtDaVSaTZfW8+ePVFQUAAvLy/ExMRYvPYdd9yB7OxsJCYmGtNqT6XSUAMHDkR8fDyWL1+O1atXo2fPnnj//feNnaUtCQ8Px549e3D33XcDMIxuPnjw4G2H1/fs2RObN29GTEyMWR+ZGgqFAv3790f//v3x/PPPIzo6Glu2bMG8efMAAD169ECPHj2wePFi9O3bF++++65ZwBMUFITIyEh8++23GDBggDE9KyvLpGZKNhzeiGYDUrYBkmux1F+EnE/N6+QHCFGzFRdLna0Gacp7rL6+Dq6grj48Bw8eNDlv3bp1wtfXV7z++uvi5MmT4ocffhD//ve/xcqVK4UQQmzdulV4eXmJ9957T5w6dUqsXr1ahISEmFz7nXfeEf7+/uLgwYPi119/FWVlZaK6ulrcddddolu3biIjI0Pk5eWJ7777Tjz77LNi7969Qggh0tPThUqlEuvXrxcnT54Uzz//vAgMDGxQH55by/HJJ58IlUolLly4IEpKSkSHDh3EoEGDxNdffy1++uknkZmZKebMmWPsy/LSSy+J0NBQsXXrVnHixAkxc+ZMERQUJMaOHWu85sCBA8VTTz1lcp+ff/5ZNG/eXIwfP17odDpx+vRp8eWXX4pHH31UVFZWiuzsbLF8+XKxd+9ecfbsWfHBBx8IpVIpPv/8c/HTTz+JRYsWiaysLHHmzBnx5ZdfipCQELFmzRohhHkfnlWrVomgoCCRnp4uTpw4IbRarfD29ha5ubl1/i2uXLli7MNjibP24WHAQ7KXnZ0tNm3aJLKzs6XOCtUjOztbvLtuncsFPEI0/j3mLgGPEIaApXv37kKpVIrg4GAxYMAA8fHHHxuPL1y4UISGhoqAgAAxceJEsWrVKpNrl5WViYceekg0a9ZMABAbNmwQQghRVFQkZs+eLSIjI4W3t7eIiooSkydPFufOnTM+d/ny5SIsLEwEBASIqVOniqSkpEYFPNXV1aJjx47iiSeeEEIIkZ+fLxITE0VYWJhQqVSiXbt2Yvr06cbvpoqKCjFr1iwRFBQkgoODhVarFX/84x/FpEmTjNe0FPAIIURubq548MEHRbNmzYSvr6+Ii4sTc+fOFdXV1SInJ0eMGDFCNG/eXKhUKhEbGyv+9re/CSGEKCgoEGPHjhUtW7YUSqVSREdHi+eff15UVVUJIcwDnqqqKrFs2TLRqlUr4e3tLbp16ya++OKLev8WrhrwKISw0KDo5IqKiqBWq1FYWHjbTm1E5EJKSoCa/gDFxYZV02WsrKwMeXl5aNu2LXx8fKTODtlZdXU17rjjDkyYMAEvvvii1Nmxm/re11J+f3O1dCIiIjs4e/Ystm3bhoEDB0Kv1+Pvf/878vLy8PDDD0udNbfkOjN7ERERuRAPDw9s3LgRd955J/r3729cN/KOO+6QOmtuiTU8REREdhAVFYXvvvtO6mzQDQx4iMh5eHgANZOwudDSEkTk/BjwEJHz8PUF9u6VOhdEJEP8CUVERESyx4CHiIiIZI8BDxE5j+vXgZgYw3b9utS5ISIZYR8eInIeQgA1ayi53pyoROTEWMNDRESys3TpUoSHh0OhUGDr1q11ppH7YMBDRESNUlBQgNmzZ6Ndu3ZQqVSIiorCAw88gJ07dzbo+ZmZmVAoFLh69WqDzi8tLcWSJUvQsWNHqFQqhIWFYfz48Th27JjJecePH8eyZcvw1ltvIT8/HyNHjrSYRu6FTVpERGS1M2fOoH///mjWrBlSU1PRtWtXVFRU4Msvv8TMmTNx4sQJm95Pr9dj6NChOHfuHFauXImEhAT88ssvWLFiBRISErBjxw5oNBoAwOnTpwEAY8aMgUKhqDON3IzDlyu1Aa6WTiRTxcUuuVp6Y7nyaukjR44UrVq1EsUWXqcrV67cdpXtmuO1t6lTp9Z5v+TkZKFQKMShQ4dM0quqqkTv3r1Fp06dRHV1tViyZInZdS2lkf0462rprOEhInJGJSV1H/P0BGqvQl3fuR4ehgkd6zvXylXpf//9d2RkZGD58uXwt/DcZs2a3baZKioqCps3b8ZDDz2EkydPIigoCL6183mLd999F8OGDUO3bt1M0j08PPD0009j8uTJOHz4MBYsWICYmBg8+uijyM/PBwAEBASYpZH7YcBDRM5DoQA6dbq5784CAuo+dt99wGef3XzcokXdw/gHDgQyM28+jokBfvvN9BwrR8SdOnUKQgjExcVZ9bzaPD09ERISAgBo0aIFmjVrVu/5ubm5GDx4sMVjNYtx5ubmonv37sZrRUREGM+xlEbuhZ2Wich5+PkBx44ZNj8/qXNDdRA3AiR79IV55513EBAQYNy++eYbSfND8sEaHiIiZ1RcXPcxT0/Tx5cu1X3urYuwnjnT6CzV6NChAxQKBY4fP46xY8fWcVvDfUWt2qOKiorbXnv06NFISEgwPm7VqhUAIDY2Fjk5ORafU9NBukOHDg3KP7kn1vAQETkjf/+6t9r9d2537q39YiydY6WQkBCMGDECb775Jkos9Am6evUqmjdvDgAmfWYOHTpkcp5SqQQAVFVVGdMCAwPRvn1741bTr2fSpEnYsWMHDh8+bHKN6upqrFq1Cp06dTLr30NUGwMeInIe168DnTsbNi4t4dTWrFmDqqoq9OnTB5s3b8aPP/6I48eP44033kDfvn3h6+sLjUaD5ORk5OTk4Ouvv8Zf//pXk2tER0dDoVDg008/xa+//oriemq1nn76afTp0wcPPPAAPvzwQ5w7dw579+7FQw89hOPHj2P9+vVs0qJ6MeAhoibT6XRIS0uDTqdr2oWEAHJyDBuXlnBqbdu2xYEDBzB48GDMnz8f8fHxGDZsGHbu3Im1a9cCAP7973+joqICvXv3xlNPPYWXXnrJ5BqtWrXCsmXLsGjRIoSHh2PWrFl13s/HxwdfffUVpk6dimeeeQbt27fHvffeCyEEMjIyEB8fb9fykutTCOF6/6sUFRVBrVajsLAQQUFBUmeHyK1ptVqkpqYaHyclJSElJaVxFyspuTk6qbi4Uc0trqSsrAx5eXlo27YtfG5tpqLbunDhAgoKCoyPIyIi0Lp1awlzRED972spv79Zw0NEjabT6UyCHQBITU1tek0P0W0UFxebBDuAYamL+prFyL0x4CGiRsvNzbUqnchW9Hq9VelEDHiIqNFiY2OtSieyFZVKZVU6EQMeImq0hIQEJCUlmaRptVqTeVSI7CEgIMBs1uSIiAgE1DdDNbk1TjxIRE2SkpKCcePGITc3F7GxsU0LdhQKIDr65j5RPVq3bo1mzZpBr9dDpVIx2KF6MeAhoiZLSEiwTa2On59NZgJ2NS44WNZp1CxBQc7DWd/PbNIiIpKIt7c3AOA6J1mUveLiYly+fNktRpGVl5cDMCwQ60xYw0NEJBFPT080a9YMl26sheXn58fZgmWooKAAv9VaoT4sLEy2q7ZXV1fj119/hZ+fH7y8nCvEcK7cEJF7Ky0FBgww7H/9tfk6UDJU88V3qb4FQMll6fV6s/mCfvvtN1y9elW2I8o8PDzQpk0bpwveGfAQkfOorgb27bu57wYUCgVatmyJFi1aNGg1cXIt//3vf6HVas3SU1JSMGbMGAlyZH9KpRIeHs7XY4YBDxHJkk6ns83IMQfx9PR0uj4P1HQxMTE4e/asxXQuJ+JYzheCERE1kVarhUajQWJiIjQajcVf2ESOwLmqnAcXDyUi52GDxUN1Oh00Go1ZenZ2Nr9kSDKuVuNoL1J+f7NJi4hkpb71vdz5i4akZbO5qqjR2KRFRLLC9b2IyBIGPETkXMLCDFsjsc8EEVnCPjxEJEvsM0HkfNiHh4jIxthngohqY5MWERERyR4DHiJyHqWlwKBBhq20VOrcEJGMsEmLiJxHdTWwe/fNfSIiG2ENDxEREckea3iIiByEI8eIpMMaHiIiB+D6XkTSYsBDRGRnOp0OqampJmmpqanQ6XQS5YjI/TDgISKys/rW9yL3otPpkJaWxmBXAgx47IBvaKIm8PMzbDLC9b0IYLOm1Bjw2Bjf0ERN4O8PlJQYNn9/qXNjM1zfi9isKT0GPDbENzQR1SUlJQXZ2dnYtGkTsrOzkZycLHWWyIHYrCk9Dku3ofre0PwlR0Rc38t9sVlTeqzhsSG+oYmaqKwMuP9+w1ZWJnVuiGyGzZrSUwghhNSZsJaUy8vfjlarNWnW0mq1rLomaqiSEiAgwLBfXCyrfjxEACeflPL7mwGPHVj7hnb3DwCREQMeIlmT8vubfXjswJp2+ltrhJKSkpCSkmKvrBGRg/EHDZFzYB8eCXFUF5G8cZoKIufBgEdCHKZIJF/8QUPkXBjwSIijuojkiz9oiJwLAx4JcZgikXzxBw2Rc+EoLSfATo1E8sRpKohMcVi6leQW8BCRfPEHDdFNUn5/W9WktXbtWnTt2hVBQUEICgpC37598cUXXwAAKioqoNVq0aVLF/j7+yMyMhKJiYm4ePGiyTX0ej1mz56NsLAw+Pv7Y/To0bhw4YLtSkRE5EQSEhIwZcoUBjtEErMq4GndujWSk5Oxb98+7Nu3D/fccw/GjBmDY8eO4fr16zhw4ACee+45HDhwAB9//DFyc3MxevRok2vMnTsXW7ZsQXp6Or799lsUFxdj1KhRqKqqsmnBiMh56XQ6pKWlmY9YKisD/vhHw8alJYjIhprcpBUSEoJXXnkF06ZNMzu2d+9e9OnTB2fPnkWbNm1QWFiI5s2bIy0tDRMnTgQAXLx4EVFRUfj8888xYsSIBt2TTVpErqveyTY50zKRrLlMk1ZtVVVVSE9PR0lJCfr27WvxnMLCQigUCjRr1gwAsH//flRUVGD48OHGcyIjIxEfH4+srKzGZoWIXATnpiEiqVi9tMSRI0fQt29flJWVISAgAFu2bEGnTp3MzisrK8OiRYvw8MMPG6O4goICKJVKBAcHm5wbHh6OgoKCOu+p1+uh1+uNj4uKiqzNNhE5gfrmpmEfFyKyJ6treDp27IhDhw4hOzsbTzzxBKZOnYqcnByTcyoqKjBp0iRUV1djzZo1t72mEAIKhaLO4ytWrIBarTZuUVFR1mabiJwA56YhIqlYHfAolUq0b98evXv3xooVK9CtWzesXr3aeLyiogITJkxAXl4etm/fbtJGFxERgfLycly5csXkmpcuXUJ4eHid91y8eDEKCwuN2/nz563NNhE5AU62SURSafJq6UIIY3NTTbDz448/YteuXQgNDTU5t1evXvD29sb27dsxYcIEAEB+fj6OHj1q1q5fm0qlgkqlampWyQ1xDhTnk5KSgnHjxvF1ISKHsirgeeaZZzBy5EhERUXh2rVrSE9PR2ZmJjIyMlBZWYnx48fjwIED+PTTT1FVVWXslxMSEgKlUgm1Wo1p06Zh/vz5CA0NRUhICBYsWIAuXbpg6NChdikgua96RwORpBISEhjoEJFDWTUsfdq0adi5cyfy8/OhVqvRtWtXaLVaDBs2DGfOnEHbtm0tPm/Xrl0YNGgQAENn5oULF+Ldd99FaWkphgwZgjVr1ljVL4fD0ul2dDodNBqNWXp2dja/aJ2ZEMD164Z9Pz+gnr59ROR6uLSElRjw0O2kpaUhMTHRLH3Tpk2YMmWKBDkiImoaOTTRu+Q8PETOjKOBiEhOtFotNBoNEhMTodFooNVqpc6Sy2HAQ7LE0UAuSq8HHnnEsNWae4vInXHCTtto8igtImfF0UAuqLISePttw/6bbwIcnUnECTtthAEPyRpHAxGRq2MTvW2wSYuIiMiJsYneNljDQ0RE5OTYRN90DHiIiIhcAJvom4YBDxER2Y0c5o4heWAfHiIisgvOHUPOhDMtE5HzEAL47TfDflgYl5ZwYVzehSzhTMtERIAhwGne3LAx2HFp9c0dQyQFBjxERGRznDuGnA0DHiJyHno9MHOmYePSEi6Nc8eQs2EfHiJyHiUlQECAYb+4GPD3lzY/1GQcpUW1Sfn9zWHpRERkN5w7hpwFm7SIiIhI9hjwEBERkewx4CEiIiLZY8BDREREsseAh4iIiGSPo7SIyHn4+gJ5eTf3iYhshAEPETkPDw8gJkbqXBCRDDHgIXIBnLyNiKhp2IeHyMlptVpoNBokJiZCo9FAq9VKnSX7KS8HFi40bOXlUueGiGSES0sQOTGdTgeNRmOWnp2dLc+aHi4tQSRrUn5/s4aHyInl5uZalU5ERJYx4CFyYrGxsValExGRZQx4iJxYQkICkpKSTNK0Wq08m7OIiOyIo7SInFxKSgrGjRvHUVrkFjgikeyFAQ+RC0hISOB//iR7Wq0WqampxsdJSUlISUmRMEckJ2zSIiIiyel0OpNgBwBSU1Oh0+kkyhHJDQMeInIevr7A0aOGjUtLuBWOSCR7Y5MWETkPDw+gc2epc0ES4IhEsjfW8BARkeQ4IpHsjTU8RA7GUSj1KC8HXn7ZsP/MM4BSKW1+yO5qfx44IpHsiUtLEDkQR6HcBpeWcCv8PLgfKb+/GfAQOYjbrYvVGAx43AY/D+6Ja2kRuQGOQiG6iZ8Hx9DpdEhLS+PwfjDgIXIYjkIhuomfB/vTarXQaDRITEyERqOBVquVOkuSYsBD5CAchUJ0Ez8P9sWJHM1xlBaRA3EUCtmSq4/44+fBfuprMnTXvzMDHiIH47pYZAtyGeHEz4N92KPJ0NUDbDZpEZHz8PEB9uwxbD4+UufGabG5gm7H1k2GcugPxGHpREQuJi0tDYmJiWbpmzZtwpQpUyTIETkrW9TK2HIKASm/v9mkRUTkYjjCiRrKFk2GcukPxCYtInIe5eXAK68YtvJyqXPjtDjCiRxJLgE2m7SIyHlwpmWruHonUnIdt3aS12q1SE5Otvo6XFrCSgx4iGSKAQ/ZGYPExrPF3459eIiIiOxMLkP5peLqUwiwDw8REckeh/ITAx4iqhMXHiS54GKlxICHiCySw0RjRDXkMtKIGo8BDxGZYfU/yQ2H8hM7LRORGckmGvPxAXbturlPZENcrNS9MeAhu+DQT9cmWfW/pycwaJB970FuzdVHGlHjsUmLbI59P1wfq/+JSG448SDZlC0XmSPpObymrqIC+Oc/DfuPPQZ4e9v/nkTkMJx4kGRDLovMkYHDq//Ly4FZswz7jzzCgIeIbIYBD9kUh35SbdXVAtcrqlBcVolifSVK9IZ/r93yuFhfieKySpQXXkPNvLfj12Zh1gPdMKhjC0nLQETywICHbKqm78eti8yxdsd1CCFQVlGNa/oKlOhvBis1Acq1G8FJsb7ixr9VKL5xruHYjefpK1FSXglrGs19y8uMAc+xi0X4pajMLmUkIvfDgIdsjkM/paGvNAQnNcGGYaswrU2xEKDcTL+5VVU7R9e+a2WVUmeBiGSCAQ/ZBYd+NkxlVfWNwKPiZg1KrQDlmjGAqTBrCrpWZqhBqQlWKqqcI0ixpRJ9ldRZICKZYMBDZKWqaoGS8ps1IxZrSW6pMbm1WaimlqWsolrq4jiUyssDgT5e8Fd5wV/phUAfLwSovBBw499gUQ6sMpz70oPx6NQ+XNoME5FsWBXwrF27FmvXrsWZM2cAAJ07d8bzzz+PkSNHAjC0/S9btgz//Oc/ceXKFSQkJODNN99E586djdfQ6/VYsGAB3nvvPZSWlmLIkCFYs2YNWrdubbtSEd1CCIHr5VXGPiiWmnHMApU69q+Xu1etg7enAv6qG4GJ6maA4q/yQuCNx/4q0+DF7Pwb6d6et5n6q6TEuPtQz9aAv7+dS0dE7sKqgKd169ZITk5G+/btAQBvv/02xowZg4MHD6Jz585ITU3Fa6+9ho0bNyI2NhYvvfQShg0bhpMnTyIwMBAAMHfuXPzvf/9Deno6QkNDMX/+fIwaNQr79++Hp6en7UtILksIAX1l9Y1mnbprT242/dxMM3l8o+nHSbqlOISHAghQeSHQxxv+Kk+zoKQmWPGvVbsSUCt4CbrxPH+VF3y8Hfi5VKmATz+9uU9EZCNNnngwJCQEr7zyCv785z8jMjISc+fONc6sq9frER4ejpSUFDz++OMoLCxE8+bNkZaWhokTJwIALl68iKioKHz++ecYMWJEg+7JiQedW3lltcnw45p+JjWje+qrZandh8WZOs86iiHg8DSpFbEUoAQaAxhDYBKo8r5Rs2LY9/H2gEKhkLo4REQmXHLiwaqqKnz44YcoKSlB3759kZeXh4KCAgwfPtx4jkqlwsCBA5GVlYXHH38c+/fvR0VFhck5kZGRiI+PR1ZWVp0Bj16vh16vNz4uKipqbLapDrU7z97sJFtlEqDU1JTUVeNSc155pXv1S/H19rRYSxLoUxO8eCPgRm1JoI9hP+BGgFKz76/yhL/SCx4eloMUw4zHRzjqjcgKXNOParM64Dly5Aj69u2LsrIyBAQEYMuWLejUqROysrIAAOHhpp0Mw8PDcfbsWQBAQUEBlEolgoODzc4pKCio854rVqzAsmXLrM2q7FUbO89WGYcfmwYrFSgpr7pRY2I6Z4ohQKkyBi+lFe7VL0Xp5WESnASoPG/UlnjdaAoydKqtHZTU7Nc+z1/pCa/b9UtpIq1WazKvUVJSElJSUup5hgurqADeecewP3kyZ1qmRnOrzw01iNUBT8eOHXHo0CFcvXoVmzdvxtSpU7F7927j8Vur0YUQt61av905ixcvxrx584yPi4qKEBUVZW3WnYIQAqUVVRY7xdZu+jEbAWQ20qfK6kndXJ2Xh8JiTUqAjxcClOYdZ0062CpN+68ovVxj3VydTmfynzYApKamYty4cfL8xVpeDjz6qGH/j39kwEON4nafG2oQqwMepVJp7LTcu3dv7N27F6tXrzb22ykoKEDLli2N51+6dMlY6xMREYHy8nJcuXLFpJbn0qVL6NevX533VKlUUEnYgbGm86xJv5RaTTpHTvyI8/m/IqBZKILCwmsdM9Sy3BqsuFO3FMWNzrO3NvXcrFmxEKBY6mDr4wWVl/v1S3HU2mSs+ie50Ol02LBhg8VjXNPPvTV5Hh4hBPR6Pdq2bYuIiAhs374dPXr0AACUl5dj9+7dxmrEXr16wdvbG9u3b8eECRMAAPn5+Th69KhZNG4LFVXV5hO51TFnSk3fFEuTvzVsUjc/4JdS4OQZm5dDCn5KT+N8KYG1aklqByYmzT+1z6vVwdZP6el2QYotOWJtMlb9k1zc+l6+Fdf0c29WBTzPPPMMRo4ciaioKFy7dg3p6enIzMxERkYGFAoF5s6di5dffhkdOnRAhw4d8PLLL8PPzw8PP/wwAECtVmPatGmYP38+QkNDERISggULFqBLly4YOnSo1Zl/4X/HUOHpU+ckb3o36zzr4+1hsRnHv3YtivJmUBJUMwHcLSOA/JVe8Kyj8yw5lr3XJmPVP8mFpfdybVzTj6wKeH755RdMmTIF+fn5UKvV6Nq1KzIyMjBs2DAAhl+GpaWlePLJJ40TD27bts04Bw8ArFq1Cl5eXpgwYYJx4sGNGzc2ag6eD/ZdgIfKz+rnORNvT4VJgFIz/0nAjdE8hqDlxtBjHy/jiB6TQOZGsHLbSd3IJdlzbTJHNZkR2Vtd7+XHH38cjz76KN/P1PR5eKRQM44/au4HkgQ8nh6GIMUbVbh4Lg+ivBTV+hJUl5dClJdiwoMPoH106xvNP7UCFGXNsOQbw5V9vKDy4mSLJB2dTgeNRmOWnp2dLc0XREkJEBBg2C8u5kzL1GBO914mi1xyHh5Xo1DA0AflRrBhEnjcmCfl5nBjC0OUa9Wm1J7U7dY2Y61Wi+TpQ6QqJpFV7N1kRuQofC/T7bh0Dc/Cd75HSHAzi6N+agc1AT5e8PP2rHNSt6biCBdydU7zHq6sBLZsMew/+CDg5Ta/ychGnOa9TBZJWcPj0gEPl5YgIiJyHWzSIiIit8QaGXIUDushIudRWQl8+KFhq6yUOjdkZ1qtFhqNBomJidBoNMYJbN2VTqdDWloadDqd1FmRJTZpOQH+wiG6gaO03AZHVZlylwlApfz+Zg2PxPgLh4jcUX1zQLmbuiYAZU2PbTHgkRDf5ETkrhyxbIqrYPDnGAx4JMQ3ORG5q5p5c2pz13lzGPw5BgMeCfFNTkTuLCUlBdnZ2di0aROys7ORnJwsdZYkweDPMdhpWWIWZ2p20w89ETstkztzhwEsnHjQSnIKeAD3eJMTNQgDHiJZ48SDbi4hIYGBDhEAKJXAhg0394lIVvbt2yfZvRnwEJHz8PYGHnlE6lwQkR3c2oXD0RjwEBGRW2J3AsexNA2Lo3GUFhE5j8pK4LPPDBuXliA74qSvjuUM062w0zIROQ92WiYH4LIWjnfr35xLSxAREdmZXCd9debFRy3NNeRoDHiIiMityHHSV1dooktJScHOnTsluz8DHiIicitym9nYldZl7N27t2T35igtIiJyOykpKRg3bpwsRmnV10TnyuWyNQY8RETkluQy6ascm+jsgU1aRERELkxuTXT2whoeInIeSiXw97/f3CeiBpFTE529cB4eIiIicggpv7/ZpEVERESyxyYtInIeVVXAN98Y9u++G/D0lDY/RCQbDHiIyHmUlQGDBxv2ubQEEdkQm7SIiIhI9hjwEBERkewx4CEiIiLZY8BDREREsseAh4iIiGSPAQ8RERHJHoelE5Hz8PYGUlNv7hMR2QgDHiJyHkolsHCh1LkgN6DT6bjulJthkxYREbkVrVYLjUaDxMREaDQaaLVaqbNEDsDFQ4nIeVRVAQcOGPZ79uTSEmRzOp0OGo3GLD07O5s1PQ7AxUOJiADD0hJ9+hi2sjKpc0MylJuba1U6yQcDHiIichuxsbFWpZN8MOAhIiK3kZCQgKSkJJM0rVbL5iw3wFFaRETkVlJSUjBu3DiO0rIXfTFQ/AtwLR+4VmDYim/8e+m8ZNliwENERG4nISGBgY61yktuBDC1ApmafWOA8wtQfq3ua+ilGyfFgIeIiMid6YtNa2FqApniX2o9Lqg/kHEBDHjI6XBCMCIiG7DUtCTDQKahGPCQU9FqtUitWVoAQFJSElJSUiTMETmUtzewZMnNfSIyZ1YjU7uvTAOblhxBFQQEhAOBETc3hRpITrr9c+2AEw+S0+CEYETk1lylRkYZCASGA4EtawUzLW8ENzfSAsIBVYDZU6X8/mYNDzmN+iYEY8BDRC6rpkbGGLy4QI1MQDgQ1LLBgYwrYMBDToMTghGqq4Hjxw37d9wBeHCqMHJidXX2vXUodnmxtPlsQo2MnDDgIadRMyFY7T48nBDMzZSWAvHxhv3iYsDfX9r8kHtyxRqZmi2gVkDjJoFMQzHgIafCCcGIyG5qBzJmNTNOVCNjKZC5tUYmMAJQ8geBNRjwkNPhhGBEZJX6Ri3VrqlxhUCGNTJ2w4CHiIick/6ahaYkFw1kWCMjOQY8RETUZFZNGKq/ZugDY9ZHJv9m+rUCoKLEMZmvCwMZWWHAQ0RETVIzYWiAEmgZ4IH5j/0Jj/9pjHkfGWerkQlqWauTbwQDGZnjxINE5DxKSoCAG/0XiouhO3qUHdidQZ1NS/kouvgj8nMPIjLQA4EqhbT5VKlvDL++dbRSeK2aGQYyUuLEg0REgGE5iQULAADPLFmCFStXGg9xmRE7qB3IWBq11ICmpSAAQWGe9s2nWSATDgRGms4tExDOQIbqxRoeInI6XGakCYQwBDK154uxNJ+MU/SRYY2Mu2ENDxFRLVxmpA631sjUNZ+MxIFMqVDCt3nMbQKZloDST9J8knthwENEzqO6Gjh3DvEBAVAAuLX6WZbLjNy2RsaZRi3VXSOTc/4KcvOL0KpjT9zZf5C0+SSygAEPETmP0lKgbVv0APDM3LlY/vrrxkMuucyIsUamjvljagIapwlkao1QMvaVaViNTKd4oJMDs0xkLfbhISLn4QqjtOqqkXHCpiWzGhmT1a8j2LREDsc+PEREFjh0mZGaQKau0Uoyq5EhcjdWBTwrVqzAxx9/jBMnTsDX1xf9+vVDSkoKOnbsaDynuLgYixYtwtatW3H58mXExMRgzpw5eOKJJ4zn6PV6LFiwAO+99x5KS0sxZMgQrFmzBq1bt7ZdyYiIgAbWyNT0kbkubV591DcDF7MamZY30xnIEFnNqoBn9+7dmDlzJu68805UVlbi2WefxfDhw5GTkwN/f8Owwaeffhq7du3Cf/7zH8TExGDbtm148sknERkZiTFjxgAA5s6di//9739IT09HaGgo5s+fj1GjRmH//v3w9LTzfA5EJA/11ciYBDKskSGiJvbh+fXXX9GiRQvs3r0bAwYMAADEx8dj4sSJeO6554zn9erVC/fddx9efPFFFBYWonnz5khLS8PEiRMBABcvXkRUVBQ+//xzjBgx4rb3ZR8eIpmq3Ycn7yBQXWh5/piagEbqGhlLgQxrZIjq5LJ9eAoLCwEAISEhxrS77roLn3zyCf785z8jMjISmZmZyM3NxerVqwEA+/fvR0VFBYYPH258TmRkJOLj45GVlWUx4NHr9dDr9cbHRUVFTck2ETmrt0ff3P/nAEAp0VIFxqalCAs1MpEMZIhcUKMDHiEE5s2bh7vuugvx8fHG9DfeeAPTp09H69at4eXlBQ8PD/zrX//CXXfdBQAoKCiAUqlEcHCwyfXCw8NRUFBg8V4rVqzAsmXLGptVInIVVWVAb2/Dvocdrm8pkGGNDJFbaHTAM2vWLPzwww/49ttvTdLfeOMNZGdn45NPPkF0dDS+/vprPPnkk2jZsiWGDh1a5/WEEFAoLP+aW7x4MebNm2d8XFRUhKioqMZmnYicVUgkcL+v9c9rSCAT2BLwbsS1iUgWGhXwzJ49G5988gm+/vprk5FVpaWleOaZZ7Blyxbcf//9AICuXbvi0KFDePXVVzF06FBERESgvLwcV65cManluXTpEvr162fxfiqVCiqVqjFZJSJXEhhh+pg1MkRkI1YFPEIIzJ49G1u2bEFmZibatm1rcryiogIVFRXw8DCti/b09ER1dTUAQwdmb29vbN++HRMmTAAA5Ofn4+jRo0hNTW1KWYjI1SU8AbQeCQS0AGI6MZAhIpuxKuCZOXMm3n33Xfz3v/9FYGCgsc+NWq2Gr68vgoKCMHDgQCxcuBC+vr6Ijo7G7t27sWnTJrz22mvGc6dNm4b58+cjNDQUISEhWLBgAbp06VJvkxcRuYGAaKBXZ8N+cTGglDY7RCQfVg1Lr6uPzYYNG/DII48AMHRKXrx4MbZt24bff/8d0dHReOyxx/D0008bn19WVoaFCxfi3XffNZl4sKH9cjgsnUimbllaAjfm9yIieZDy+5traRGR82DAQyRrLjsPDxERWU+n0znfoqhEMmePmS6IiKgOWq0WGo0GiYmJ0Gg00Gq1UmeJyC0w4CEichCdTmc2GjU1NRU6nU6iHBG5DwY8REQOkpuba1U6EdkO+/AQkfPw8gKmTr25LzOxsbFWpROR7bCGh4ich0oFbNxo2GQ4u3pCQgKSkpJM0rRaLTsuEzkAh6UTETkYR2mRu+KwdCIiABACuH7dsO/nB9Qx2amrS0hIYKBD5GBs0iIi53H9umHiwYCAm4EPEZENMOAhIiIi2WPAQ0RERLLHgIeIiIhkjwEPERERyR4DHiIiIpI9BjxEREQke5yHh4ich6cnMH78zX0iIhthwENEzsPHB/jwQ6lzQUQyxCYtIiIikj0GPERERCR7DHiIyHmUlBjWz1IoDPtERDbCgIeIiIhkjwEPERERyR4DHiIiIpI9BjxEREQkewx4iIiISPYY8BAREZHscaZlInIenp7Afffd3CcishEGPETkPHx8gM8+kzoXRCRDbNIiIiIi2WMNDxERNZlOp0Nubi5iY2ORkJAgdXaIzLCGh4icR0kJ4O9v2Li0hMvQarXQaDRITEyERqOBVquVOktEZhjwEJFzuX7dsJFL0Ol0SE1NNUlLTU2FTqeTKEdEljHgoTrpdDqkpaXxPy4iqlNubq5V6URSYcBDFrGKmogaIjY21qp0Iqkw4CEzrKImooZKSEhAUlKSSZpWq2XHZXI6HKVFZuqrouZ/YkR0q5SUFIwbN46jtMipMeAhM6yiJiJrJSQkMNAhp8YmLTLDKmqSjIcHMHCgYfPgf09EZDsKIYSQOhPWKioqglqtRmFhIYKCgqTOjmxxIjEiIrIlKb+/2aRFdWIVNRERyQXrjImIiEj2GPAQkfMoKQGaNzdsXFqCiGyITVpE5Fx++03qHBCRDLGGh4iIiGSPAQ8RERHJHpu0iIiIqMFcdcoS1vAQERFRg7jywtIMeIiIiOi2XH1haQY8ROQ8PDyA3r0NG5eWIHIq9S0s7QrYh4eInIevL7B3r9S5ICILXH1haf6EIiIiotty9YWluXgoERERNVhTRmlx8VAiIgC4fh3o1Mmwn5MD+PlJmx8iMuOqC0sz4CEi5yEEcPbszX0iIhthwENETs9VJzojIufBTstE5NRceaIzInIeDHiIyGm5+kRnROQ8GPAQkdNy9YnOiMh5sA8PETktqSc6Y98hIvlgDQ8ROQ+FwjAsvVMnQKGQdKIz9h0ikhdOPEhNxl/BZG+Ofo/pdDpoNBqz9OzsbL7HiZqAEw+Sy9JqtSadSpOSkpCSkiJhjkiOHD3RWX19hxjwELkmNmlRo3EEDcmV1H2HiMj2rAp4VqxYgTvvvBOBgYFo0aIFxo4di5MnT5qdd/z4cYwePRpqtRqBgYHQaDQ4d+6c8bher8fs2bMRFhYGf39/jB49GhcuXGh6acihOIKGbO76daBzZ8N2/bpk2XD1RRKJyJxVAc/u3bsxc+ZMZGdnY/v27aisrMTw4cNRUlJiPOf06dO46667EBcXh8zMTBw+fBjPPfccfHx8jOfMnTsXW7ZsQXp6Or799lsUFxdj1KhRqKqqsl3JyO74K5hsTgjDGlo5OZIvLZGSkoLs7Gxs2rQJ2dnZSE5OljQ/Op0OaWlprEElaizRBJcuXRIAxO7du41pEydOFP/3f/9X53OuXr0qvL29RXp6ujHt559/Fh4eHiIjI6NB9y0sLBQARGFhYeMzTzaRlJQkABg3rVYrdZbIlRUXC2EIdQz7JIQw/5wlJSVJnSWiRpHy+7tJfXgKCwsBACEhIQCA6upqfPbZZ4iNjcWIESPQokULJCQkYOvWrcbn7N+/HxUVFRg+fLgxLTIyEvHx8cjKyrJ4H71ej6KiIpONnIOz/Qomkhv2lSOyjUYHPEIIzJs3D3fddRfi4+MBAJcuXUJxcTGSk5Nx7733Ytu2bXjwwQcxbtw47N69GwBQUFAApVKJ4OBgk+uFh4ejoKDA4r1WrFgBtVpt3KKiohqbbbKDhIQETJkyhf0biOzAGfrKsTmN5KDRAc+sWbPwww8/4L333jOmVVdXAwDGjBmDp59+Gt27d8eiRYswatQo/OMf/6j3ekIIKBQKi8cWL16MwsJC43b+/PnGZpuIyKVI3VeOEzCSXDQq4Jk9ezY++eQT7Nq1C61btzamh4WFwcvLC506dTI5/4477jCO0oqIiEB5eTmuXLlics6lS5cQHh5u8X4qlQpBQUEmGxGRO5ByxBib00hOrAp4hBCYNWsWPv74Y3z11Vdo27atyXGlUok777zTbKh6bm4uoqOjAQC9evWCt7c3tm/fbjyen5+Po0ePol+/fo0tR5OxypbICSgUQHS0YaujxtcdSdVXzhma04hsxaqZlmfOnIl3330X//3vfxEYGGjsc6NWq+Hr6wsAWLhwISZOnIgBAwZg8ODByMjIwP/+9z9kZmYaz502bRrmz5+P0NBQhISEYMGCBejSpQuGDh1q29I1EGcLJnISfn7AmTNS58IpOXq2aUD65jQim7JmSBdqDYusvW3YsMHkvPXr14v27dsLHx8f0a1bN7F161aT46WlpWLWrFkiJCRE+Pr6ilGjRolz5841OB+1h7VlZ2eLTZs2iezsbGuKYpSdnW2xTI29HhGRnHDqCbIlKYelu/TioU899RRWr15tTG9MzUxaWhoSExPN0jdt2oQpU6Y0Oa9ERK6OCwSTrUi5eKhLBzyWWLuaMVdFJnIipaXAgAGG/a+/Bm40lRORPEgZ8Mhu8VBrO9NxzRwiJ1JdDezbZ9huTHNBRGQLVnVadgWN6UyXkpKCcePGscqWiIhIply6huepp54yedyUmhnOFkxERCRfLt2Hp7CwEMePH2fNDJFclJQAAQGG/eJiwN9f2vwQkU1J2YfH5Zu0pJibgojIXXCEFsmFSzdpERGR/XAdLZITBjxE5FzCwgwbSYrraJHcMOAhIufh7w/8+qthY/8dSXEdLZIbBjxERGSG62iR3DDgISIiM5yUleTG5YelO3pYGxHZUWkpMHKkYf+LL7i0hBPgKC2yJa6lZSUGPEQyxXl4iGSNa2kRERER2REDHiIiIpI9BjxEREQkewx4iIiISPYY8BAREZHsufziobbGIZhEEvPzkzoHRLLA7zNTrOGphQvlEUnM398wNL2khEPSiZqA32fmOA/PDTqdDhqNxiw9OzubkTEREbkMZ/4+4zw8ToAL5RERkRzw+8wyBjw3cKE8IidQVgbcf79hKyuTOjdELonfZ5Yx4LmBC+UROYGqKuDzzw1bVZXUuSFySfw+s4x9eG7BXu1EEuJaWkQ244zfZ1w81EpcPJRIphjwEMkaOy0TERER2REDHiIiIpI9BjxEREQkey65tERNt6OioiKJc0JENlVScnO/qIgjtYhkpuZ7W4ruwy4Z8Fy+fBkAEBUVJXFOiMhuIiOlzgER2cnly5ehVqsdek+XDHhCQkIAAOfOnXP4H0xKRUVFiIqKwvnz591qdBrLzXK7A5ab5XYHhYWFaNOmjfF73JFcMuDx8DB0PVKr1W71RqkRFBTEcrsRltu9sNzuxV3LXfM97tB7OvyORERERA7GgIeIiIhkzyUDHpVKhSVLlkClUkmdFYdiuVlud8Bys9zugOV2fLldcmkJIiIiImu4ZA0PERERkTUY8BAREZHsMeAhIiIi2WPAQ0RERLLnFAHPlStXMGXKFKjVaqjVakyZMgVXr16t9zkff/wxRowYgbCwMCgUChw6dMjsHL1ej9mzZyMsLAz+/v4YPXo0Lly40OR720pj7i2EwNKlSxEZGQlfX18MGjQIx44dMx4/c+YMFAqFxe3DDz80nhcTE2N2fNGiRfYqqgl7lBsABg0aZFamSZMmNfnetmKPcv/++++YPXs2OnbsCD8/P7Rp0wZz5sxBYWGhyXUc+XqvWbMGbdu2hY+PD3r16oVvvvmm3vN3796NXr16wcfHB+3atcM//vEPs3M2b96MTp06QaVSoVOnTtiyZUuT72trti73unXrcPfddyM4OBjBwcEYOnQo9uzZY3LO0qVLzV7XiIgIm5etPrYu98aNGy3+/1VWVtak+9qarctt6f8vhUKB+++/33iOq73e+fn5ePjhh9GxY0d4eHhg7ty5Fs9z2OdbOIF7771XxMfHi6ysLJGVlSXi4+PFqFGj6n3Opk2bxLJly8S6desEAHHw4EGzc2bMmCFatWoltm/fLg4cOCAGDx4sunXrJiorK5t0b1tpzL2Tk5NFYGCg2Lx5szhy5IiYOHGiaNmypSgqKhJCCFFZWSny8/NNtmXLlgl/f39x7do143Wio6PFCy+8YHJe7eP2ZI9yCyHEwIEDxfTp003KdPXq1Sbf21bsUe4jR46IcePGiU8++UScOnVK7Ny5U3To0EE89NBDJtdx1Oudnp4uvL29xbp160ROTo546qmnhL+/vzh79qzF83/66Sfh5+cnnnrqKZGTkyPWrVsnvL29xUcffWQ8JysrS3h6eoqXX35ZHD9+XLz88svCy8tLZGdnN/q+tmaPcj/88MPizTffFAcPHhTHjx8Xjz76qFCr1eLChQvGc5YsWSI6d+5s8rpeunTJ7uWtYY9yb9iwQQQFBZn9P9aU+9qaPcp9+fJlk/IePXpUeHp6ig0bNhjPcbXXOy8vT8yZM0e8/fbbonv37uKpp54yO8eRn2/JA56cnBwBwKRw33//vQAgTpw4cdvn5+XlWQx4rl69Kry9vUV6erox7eeffxYeHh4iIyPDJvduisbcu7q6WkRERIjk5GRjWllZmVCr1eIf//hHnffq3r27+POf/2ySFh0dLVatWtW0QjSCPcs9cOBAix+optzbVhz5en/wwQdCqVSKiooKY5qjXu8+ffqIGTNmmKTFxcWJRYsWWTw/KSlJxMXFmaQ9/vjjQqPRGB9PmDBB3HvvvSbnjBgxQkyaNKnR97U1e5T7VpWVlSIwMFC8/fbbxrQlS5aIbt26NT7jTWSPcm/YsEGo1Wqb3tfWHPF6r1q1SgQGBori4mJjmqu93rXV9f+zIz/fkjdpff/991Cr1UhISDCmaTQaqNVqZGVlNfq6+/fvR0VFBYYPH25Mi4yMRHx8vPG69rp3QzTm3nl5eSgoKDApk0qlwsCBA+t8zv79+3Ho0CFMmzbN7FhKSgpCQ0PRvXt3LF++HOXl5U0s1e3Zu9zvvPMOwsLC0LlzZyxYsADXrl1r0r1txVGvN2BYnC8oKAheXqZL5dn79S4vL8f+/ftN8gsAw4cPrzO/33//vdn5I0aMwL59+1BRUVHvOTXXbMx9bcle5b7V9evXUVFRYbbo4o8//ojIyEi0bdsWkyZNwk8//dSE0jScPctdXFyM6OhotG7dGqNGjcLBgwebdF9bctTrvX79ekyaNAn+/v4m6a70ejeEIz/fki8eWlBQgBYtWpilt2jRAgUFBU26rlKpRHBwsEl6eHi48br2undD82ftvWvSw8PDTdLDw8Nx9uxZi89Zv3497rjjDvTr188k/amnnkLPnj0RHByMPXv2YPHixcjLy8O//vWvxhSnwexZ7smTJ6Nt27aIiIjA0aNHsXjxYhw+fBjbt29v9L1txVGv9+XLl/Hiiy/i8ccfN0l3xOv922+/oaqqymJ+6yujpfMrKyvx22+/oWXLlnWeU3PNxtzXluxV7lstWrQIrVq1wtChQ41pCQkJ2LRpE2JjY/HLL7/gpZdeQr9+/XDs2DGEhobaoHR1s1e54+LisHHjRnTp0gVFRUVYvXo1+vfvj8OHD6NDhw5u8Xrv2bMHR48exfr1603SXe31bghHfr7tFvAsXboUy5Ytq/ecvXv3AgAUCoXZMSGExfSmuvW6tr63I8p96/G6nlNaWop3330Xzz33nNmxp59+2rjftWtXBAcHY/z48cZaAGs5Q7mnT59u3I+Pj0eHDh3Qu3dvHDhwAD179mzSveviDOWuUVRUhPvvvx+dOnXCkiVLTI7Z+vW2RX7rO//W9IZc09r72po9yl0jNTUV7733HjIzM+Hj42NMHzlypHG/S5cu6Nu3L/7whz/g7bffxrx58xpVDmvZutwajQYajcZ4vH///ujZsyf+9re/4Y033mj0fW3Nnq/3+vXrER8fjz59+piku+Lrbatr2uK+dgt4Zs2aZTZC5lYxMTH44Ycf8Msvv5gd+/XXX80iOmtERESgvLwcV65cManluXTpkrG2IyIiwub3tme5a3rjFxQUmPwiuHTpksXnfPTRR7h+/ToSExNvm++a/2BOnTrVqC9AZyp3jZ49e8Lb2xs//vgjevbsKevX+9q1a7j33nsREBCALVu2wNvbu948NfX1tiQsLAyenp5mv7rqe50iIiIsnu/l5WXMV13n1FyzMfe1JXuVu8arr76Kl19+GTt27EDXrl3rzYu/vz+6dOmCH3/8sRElsY69y13Dw8MDd955p7FMcn+9r1+/jvT0dLzwwgu3zYuzv94N4cjPt9368ISFhSEuLq7ezcfHB3379kVhYaHJcEudTofCwkKzZhhr9OrVC97e3sbmDMAwRO7o0aPG69rj3vYsd01zTe0ylZeXY/fu3Rafs379eowePRrNmze/bb5r2sgtVaU3hDOVu8axY8dQUVFhLJNcX++ioiIMHz4cSqUSn3zyiUkNQF2a+npbolQq0atXL5P8AsD27dvrLGPfvn3Nzt+2bRt69+5tDNrqOqfmmo25ry3Zq9wA8Morr+DFF19ERkYGevfufdu86PV6HD9+3Kava13sWe7ahBA4dOiQsUxyfr0B4IMPPoBer8f//d//3TYvzv56N4RDP99WdXG2k3vvvVd07dpVfP/99+L7778XXbp0MRuu27FjR/Hxxx8bH1++fFkcPHhQfPbZZwKASE9PFwcPHjQZvjhjxgzRunVrsWPHDnHgwAFxzz33WByWfrt720tjyp2cnCzUarX4+OOPxZEjR8Sf/vQns+HZQgjx448/CoVCIb744guz+2ZlZYnXXntNHDx4UPz000/i/fffF5GRkWL06NH2Kegt7FHuU6dOiWXLlom9e/eKvLw88dlnn4m4uDjRo0cPWb/eRUVFIiEhQXTp0kWcOnXKZLhqTbkd+XrXDB9dv369yMnJEXPnzhX+/v7izJkzQgghFi1aJKZMmWI8v2a47tNPPy1ycnLE+vXrzYbrfvfdd8LT01MkJyeL48ePi+Tk5DqHrdZ1X3uzR7lTUlKEUqkUH330UZ3TCcyfP19kZmaKn376SWRnZ4tRo0aJwMBAly730qVLRUZGhjh9+rQ4ePCgePTRR4WXl5fQ6XQNvq8rlrvGXXfdJSZOnGjxvq72egshxMGDB8XBgwdFr169xMMPPywOHjwojh07ZjzuyM+3UwQ8ly9fFpMnTxaBgYEiMDBQTJ48WVy5csXkHAAm8xFs2LBBADDblixZYjyntLRUzJo1S4SEhAhfX18xatQoce7cOavvbS+NKXd1dbVYsmSJiIiIECqVSgwYMEAcOXLE7NqLFy8WrVu3FlVVVWbH9u/fLxISEoRarRY+Pj6iY8eOYsmSJaKkpMTWRbTIHuU+d+6cGDBggAgJCRFKpVL84Q9/EHPmzBGXL1+2+t72Yo9y79q1y+LnAIDIy8sTQjj+9X7zzTdFdHS0UCqVomfPnmL37t3GY1OnThUDBw40OT8zM1P06NFDKJVKERMTI9auXWt2zQ8//FB07NhReHt7i7i4OLF582ar7usIti53dHT0bf+Pq5mXydvbW0RGRopx48aZfJk4gq3LPXfuXNGmTRuhVCpF8+bNxfDhw0VWVpZV93UEe7zPT548KQCIbdu2WbynK77elt7D0dHRJuc46vOtuJEhIiIiItmSfB4eIiIiIntjwENERESyx4CHiIiIZI8BDxEREckeAx4iIiKSPQY8REREJHsMeIiIiEj2GPAQERGR7DHgISIiItljwENERESyx4CHiIiIZI8BDxEREcne/wODLlBQ1K5vtQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(df.pov_index, df.conspcfood_t1, s=10, color='black')  # s controls the size\n",
    "neg_xval = np.linspace(-0.1, 0, 100)\n",
    "neg_line = model_control.predict(neg_xval.reshape(-1, 1))\n",
    "pos_xval = np.linspace(0, 0.1, 100)\n",
    "pos_line = model_treated.predict(pos_xval.reshape(-1, 1))\n",
    "plt.plot(neg_xval, neg_line, linewidth=3, label=\"Control Regression\")\n",
    "plt.plot(pos_xval, pos_line, linewidth=3, label=\"Treated Regression\")\n",
    "plt.axvline(x=0, color='red', linestyle='--', label=\"Cut-Off\")\n",
    "plt.legend()\n",
    "plt.xlim(-0.1, 0.1)\n",
    "plt.ylim(250, 350)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "W6MN7zBmAxn9"
   },
   "source": [
    "We can repeat the estimation using the `rdd` package, which yields us an estimate as well as a confidence band calculated according to the formulas presented in the book. We look at all four targets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "Xkqz2Bk_56v4"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "289"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.pov_index.duplicated().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "id": "6rLo9c_YGWIq"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Food T_1</th>\n",
       "      <td>-22.161738</td>\n",
       "      <td>27.456747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Non-Food T_1</th>\n",
       "      <td>-9.135265</td>\n",
       "      <td>21.932340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Food T_2</th>\n",
       "      <td>54.956943</td>\n",
       "      <td>48.120411</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Non-Food T_2</th>\n",
       "      <td>43.806595</td>\n",
       "      <td>32.352294</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   LATE       s.e.\n",
       "Food T_1     -22.161738  27.456747\n",
       "Non-Food T_1  -9.135265  21.932340\n",
       "Food T_2      54.956943  48.120411\n",
       "Non-Food T_2  43.806595  32.352294"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = []\n",
    "for outcome in [\"conspcfood_t1\", \"conspcnonfood_t1\", \"conspcfood_t2\", \"conspcnonfood_t2\"]:\n",
    "    rdd_result = rdrobust(x=df.pov_index, y=df[outcome], rho=1, masspoints=\"off\")\n",
    "    result.append([rdd_result.coef.iloc[0].values[0], rdd_result.se.iloc[2].values[0]])\n",
    "res_dataframe = pd.DataFrame(result, columns=[\"LATE\", \"s.e.\"],\n",
    "                             index=[\"Food T_1\", \"Non-Food T_1\", \"Food T_2\", \"Non-Food T_2\"])\n",
    "res_dataframe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BzzCc3oWZycJ"
   },
   "source": [
    "While the effects in the first year after the intervention are negative, we observe significant positive effects in the second year after an individual or household was accepted in the *Progresa* program. This is in accordance to the previous analysis of this dataset. One possible explanation for this is that the program households have more money and can thus afford more. This was the desired effect of the program to combat hunger and extreme poverty."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hDEf53bE-Aki"
   },
   "source": [
    "## Estimation with Covariates\n",
    "\n",
    "For the identification and estimation of the average treatment effect at the cutoff value no covariate information is required except the running variable, but nevertheless in many applications additional covariates are collected which might be exploited for the analysis to improve the efficiency of the estimates.\n",
    "\n",
    "The standard approach is simply to take up the regressors in the weighted least squares regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "JRdUQ8gcsGCg"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-23.781010315693663"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_treated, model_control = LinearRegression(), LinearRegression()\n",
    "model_treated.fit(y=df.loc[df.pov_index > 0, \"conspcfood_t1\"].values.reshape(-1, 1),\n",
    "                  X=df.loc[df.pov_index > 0, [\"pov_index\", \"hhownhouse\", \"headage\", \"heademp\", \"headeduc\"]],\n",
    "                  sample_weight=weights[df.pov_index > 0])\n",
    "model_control.fit(y=df.loc[df.pov_index < 0, \"conspcfood_t1\"].values.reshape(-1, 1),\n",
    "                  X=df.loc[df.pov_index < 0, [\"pov_index\", \"hhownhouse\", \"headage\", \"heademp\", \"headeduc\"]],\n",
    "                  sample_weight=weights[df.pov_index < 0])\n",
    "pred_t = model_treated.predict(pd.DataFrame({\"pov_index\": cutoff,\n",
    "                                             \"hhownhouse\": np.average(df.loc[df.pov_index > 0, \"hhownhouse\"],\n",
    "                                                                      weights=weights[df.pov_index > 0]),\n",
    "                                             \"headage\": np.average(df.loc[df.pov_index > 0, \"headage\"],\n",
    "                                                                   weights=weights[df.pov_index > 0]),\n",
    "                                             \"heademp\": np.average(df.loc[df.pov_index > 0, \"heademp\"],\n",
    "                                                                   weights=weights[df.pov_index > 0]),\n",
    "                                             \"headeduc\": np.average(df.loc[df.pov_index > 0, \"headeduc\"],\n",
    "                                                                    weights=weights[df.pov_index > 0])},\n",
    "                                            index=[0]))\n",
    "pred_c = model_control.predict(pd.DataFrame({\"pov_index\": cutoff,\n",
    "                                             \"hhownhouse\": np.average(df.loc[df.pov_index < 0, \"hhownhouse\"],\n",
    "                                                                      weights=weights[df.pov_index < 0]),\n",
    "                                             \"headage\": np.average(df.loc[df.pov_index < 0, \"headage\"],\n",
    "                                                                   weights=weights[df.pov_index < 0]),\n",
    "                                             \"heademp\": np.average(df.loc[df.pov_index < 0, \"heademp\"],\n",
    "                                                                   weights=weights[df.pov_index < 0]),\n",
    "                                             \"headeduc\": np.average(df.loc[df.pov_index < 0, \"headeduc\"],\n",
    "                                                                    weights=weights[df.pov_index < 0])},\n",
    "                                            index=[0]))\n",
    "treatment_effect = pred_t - pred_c\n",
    "treatment_effect[0][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Mlc6ZqDuMfl8"
   },
   "source": [
    "Including these selected covariates does not have a significant impact on the LATE estimation. Again, we can also use `rdd` to repeat the estimation with all other outcomes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "GZdzDHEdOt-j"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Food T_1</th>\n",
       "      <td>-22.161738</td>\n",
       "      <td>27.456747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Non-Food T_1</th>\n",
       "      <td>-9.135265</td>\n",
       "      <td>21.932340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Food T_2</th>\n",
       "      <td>54.956943</td>\n",
       "      <td>48.120411</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Non-Food T_2</th>\n",
       "      <td>43.806595</td>\n",
       "      <td>32.352294</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   LATE       s.e.\n",
       "Food T_1     -22.161738  27.456747\n",
       "Non-Food T_1  -9.135265  21.932340\n",
       "Food T_2      54.956943  48.120411\n",
       "Non-Food T_2  43.806595  32.352294"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "id": "OYQuZcvjyYx6"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "      <th>% reduction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Food T_1</th>\n",
       "      <td>-27.137309</td>\n",
       "      <td>21.871141</td>\n",
       "      <td>-20.343291</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Non-Food T_1</th>\n",
       "      <td>-5.068397</td>\n",
       "      <td>20.534658</td>\n",
       "      <td>-6.372699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Food T_2</th>\n",
       "      <td>57.286722</td>\n",
       "      <td>44.369975</td>\n",
       "      <td>-7.793856</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Non-Food T_2</th>\n",
       "      <td>47.665945</td>\n",
       "      <td>29.601985</td>\n",
       "      <td>-8.501124</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   LATE       s.e.  % reduction\n",
       "Food T_1     -27.137309  21.871141   -20.343291\n",
       "Non-Food T_1  -5.068397  20.534658    -6.372699\n",
       "Food T_2      57.286722  44.369975    -7.793856\n",
       "Non-Food T_2  47.665945  29.601985    -8.501124"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = []\n",
    "for outcome in [\"conspcfood_t1\", \"conspcnonfood_t1\", \"conspcfood_t2\", \"conspcnonfood_t2\"]:\n",
    "    rdd_result = rdrobust(x=df.pov_index, y=df[outcome], rho=1, masspoints=\"off\",\n",
    "                          covs=df.iloc[:, [0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 18, 21]])\n",
    "    result.append([rdd_result.coef.iloc[0].values[0], rdd_result.se.iloc[2].values[0]])\n",
    "res_dataframe_adj = pd.DataFrame(result, columns=[\"LATE\", \"s.e.\"],\n",
    "                                 index=[\"Food T_1\", \"Non-Food T_1\", \"Food T_2\", \"Non-Food T_2\"])\n",
    "res_dataframe_adj[\"% reduction\"] = (res_dataframe_adj[\"s.e.\"] - res_dataframe[\"s.e.\"]) * 100 / res_dataframe[\"s.e.\"]\n",
    "res_dataframe_adj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5q8S0wNhabWy"
   },
   "source": [
    "Overall, the adjustment by only a few covariates has not changed the estimated coefficient much from the result without covariates. However, including covariates does reduce the standard deviation of the estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9U8UkHmv-D-0"
   },
   "source": [
    "## Estimation using ML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NiYSglH9E0Er"
   },
   "source": [
    "As discussed in the book, including many covariates in RDD estimation can be beneficial for multiple reasons:\n",
    "1. **Efficiency and power improvements**: As in randomized control trials, using covariates can increase efficiency and improve power.\n",
    "2. **Auxiliary information**: In RDD the score determines the treatment assignment and measurement errors in the running variable can distort the results. Additional covariates can be exploited to overcome these issues or to deal with missing data problems.\n",
    "3. **Treatment effect heterogeneity**: Covariates can be used to define subgroups in which the treatment effects differ.\n",
    "4. **Other parameters of interest and extrapolation**: As the identified treatment effect in RDD is local at the cutoff, additional covariates might help for extrapolation of the treatment effects or identify other causal parameters.\n",
    "\n",
    "However, including a high number of covariates also comes with additional challenges, such as variables selection, non-linearities or interactions between covariates. The best way to overcome these is the use of modern ML methods.\n",
    "\n",
    "There are multiple ways to implement the estimators presented in the book, we will closely follow the analysis of [Noack et al. (2023)](https://arxiv.org/pdf/2107.07942.pdf). We set up running variable and outcome as above. The baseline covariates will be all the other variables in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "n2uoMwzkCq4P"
   },
   "outputs": [],
   "source": [
    "# Running Variable and Outcome\n",
    "investigated_outcome = \"conspcfood_t1\"\n",
    "df_ml = df.rename(columns={\"pov_index\": \"X\", investigated_outcome: \"Y\"})\n",
    "\n",
    "# Baseline covariates including consumption\n",
    "b_covs = df_ml.columns[[0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 18, 21]]\n",
    "\n",
    "# Fixed effects for localities\n",
    "i_fe = pd.get_dummies(df_ml['clus'], drop_first=True)\n",
    "\n",
    "# Flexible covariates including localities indicators\n",
    "f_covs = patsy.dmatrix('~ (' + ' + '.join(b_covs) + ')**2', data=df_ml, return_type='dataframe')\n",
    "\n",
    "# Dropping the intercept column that is automatically added by patsy\n",
    "f_covs = f_covs.iloc[:, 1:]\n",
    "\n",
    "Z_lasso = pd.concat([i_fe, f_covs], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2n8yvua4Ns_A"
   },
   "source": [
    "We will use the package `rdrobust` for the RD estimation. Before starting the DML procedure, we have to estimate a bandwidth to restrict the samples in the first stage estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "id": "VIO-PQEtOKob"
   },
   "outputs": [],
   "source": [
    "h_fs = rdrobust(y=df_ml.Y, x=df_ml.X, masspoints=\"off\").bws.values[1, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nm8BC6JTQnV7"
   },
   "source": [
    "The next chunk sets up the crossfitting and estimates the function $\\eta(Z)$, which we will use to adjust $Y$ for the second stage. We use Random Forest, a Boosting implementation, Linear Regression and Lasso with both a baseline and flexible covariate structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "id": "y-tGMe5iQhVd"
   },
   "outputs": [],
   "source": [
    "def first_stage():\n",
    "    # Set up the cross-fitting\n",
    "    n = df_ml.shape[0]\n",
    "    Kf = 5  # Number of folds\n",
    "    sampleframe = np.repeat(range(0, Kf), math.ceil(n / Kf))\n",
    "    cfgroup = np.random.choice(sampleframe, size=n, replace=False)\n",
    "\n",
    "    # Matrix to store eta predictions\n",
    "    eta_fit = np.empty((n, 5))\n",
    "\n",
    "    # Create vector of observations to be considered in the first stage model\n",
    "    weights = np.abs(df_ml.X) < h_fs\n",
    "\n",
    "    for k in range(Kf):\n",
    "        fold = (cfgroup == k)\n",
    "\n",
    "        data_treated = df_ml.loc[(df_ml.X > 0) & ~fold & (weights > 0)]\n",
    "        data_control = df_ml.loc[(df_ml.X < 0) & ~fold & (weights > 0)]\n",
    "\n",
    "        data_fold = df_ml.loc[fold]\n",
    "\n",
    "        rf1 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        rf0 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 0] = (rf1.predict(data_fold[b_covs]) + rf0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lgbm1 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lgbm0 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 1] = (lgbm1.predict(data_fold[b_covs]) + lgbm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lm1 = LinearRegression()\n",
    "        lm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lm0 = LinearRegression()\n",
    "        lm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 2] = (lm1.predict(data_fold[b_covs]) + lm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        las_base1 = LassoCV()\n",
    "        las_base1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        las_base0 = LassoCV()\n",
    "        las_base0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 3] = (las_base1.predict(data_fold[b_covs]) + las_base0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        X_flex_treated = pd.concat([Z_lasso.loc[data_treated.index], data_treated[b_covs]], axis=1)\n",
    "        X_flex_control = pd.concat([Z_lasso.loc[data_control.index], data_control[b_covs]], axis=1)\n",
    "        X_flex_fold = pd.concat([Z_lasso.loc[data_fold.index], data_fold[b_covs]], axis=1)\n",
    "        X_flex_treated.columns = X_flex_treated.columns.astype(str)\n",
    "        X_flex_control.columns = X_flex_control.columns.astype(str)\n",
    "        X_flex_fold.columns = X_flex_fold.columns.astype(str)\n",
    "        las_flex1 = LassoCV()\n",
    "        las_flex1.fit(y=data_treated.Y, X=X_flex_treated)\n",
    "        las_flex0 = LassoCV()\n",
    "        las_flex0.fit(y=data_control.Y, X=X_flex_control)\n",
    "        eta_fit[fold, 4] = (las_flex1.predict(X_flex_fold) + las_flex0.predict(X_flex_fold)) / 2\n",
    "\n",
    "    return eta_fit\n",
    "\n",
    "\n",
    "eta_fit = first_stage()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ybTRUohWi_xE"
   },
   "source": [
    "With the estimated $\\hat{\\eta}(Z)$ we can correct for confounding in $Y$ and now run the RD estimation as second stage again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "id": "WdJkfePmx4iN"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>-17.290232</td>\n",
       "      <td>20.260055</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient Boosting</th>\n",
       "      <td>-18.029527</td>\n",
       "      <td>21.068436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Regression</th>\n",
       "      <td>-24.438212</td>\n",
       "      <td>20.175355</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Baseline</th>\n",
       "      <td>-18.652900</td>\n",
       "      <td>19.925989</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Flexible</th>\n",
       "      <td>-20.075361</td>\n",
       "      <td>19.858670</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        LATE       s.e.\n",
       "Random Forest     -17.290232  20.260055\n",
       "Gradient Boosting -18.029527  21.068436\n",
       "Linear Regression -24.438212  20.175355\n",
       "Lasso Baseline    -18.652900  19.925989\n",
       "Lasso Flexible    -20.075361  19.858670"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "methods = [\"Random Forest\", \"Gradient Boosting\", \"Linear Regression\",\n",
    "           \"Lasso Baseline\", \"Lasso Flexible\"]\n",
    "\n",
    "\n",
    "def second_stage(eta_fit):\n",
    "    adj_results = []\n",
    "    for i in range(len(methods)):\n",
    "        M_Y = df_ml.Y - eta_fit[:, i]\n",
    "        rd_call = rdrobust(y=M_Y, x=df_ml.X, masspoints=\"off\")\n",
    "        adj_results.append([rd_call.coef.iloc[0].values[0],\n",
    "                            rd_call.se.iloc[2].values[0]])\n",
    "    return adj_results\n",
    "\n",
    "\n",
    "adj_frame = pd.DataFrame(second_stage(eta_fit), columns=[\"LATE\", \"s.e.\"],\n",
    "                         index=methods)\n",
    "adj_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dWRrjcieIKwk"
   },
   "source": [
    "Finally, we create a small simulation study with only $R=20$ repetitions to show the variance reducing effect of the inclusion of ML-based estimators for the covariates. The next block runs up to ten minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "id": "3DxmPT9LIJ3E"
   },
   "outputs": [],
   "source": [
    "estimates = [adj_frame[\"LATE\"].values]\n",
    "std_err = [adj_frame[\"s.e.\"].values]\n",
    "R = 19\n",
    "\n",
    "for i in range(R):\n",
    "    eta_fit = first_stage()\n",
    "    adj_results = np.array(second_stage(eta_fit))\n",
    "    estimates.append(adj_results[:, 0])\n",
    "    std_err.append(adj_results[:, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "W_ybpwiMS1MH"
   },
   "source": [
    "We aggregate the median of the estimates, the mean of the standard errors and also calculate the mean reduction of standard error compared to the \"no covariates\" estimation. We see, that including covariates can reduce the standard error of estimation around 15-20%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "id": "bzp90cR4I6RL"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "      <th>% reduction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>-19.591159</td>\n",
       "      <td>20.252981</td>\n",
       "      <td>-26.236779</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient Boosting</th>\n",
       "      <td>-21.473448</td>\n",
       "      <td>21.536645</td>\n",
       "      <td>-21.561559</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Regression</th>\n",
       "      <td>-27.272582</td>\n",
       "      <td>20.229335</td>\n",
       "      <td>-26.322902</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Baseline</th>\n",
       "      <td>-22.005207</td>\n",
       "      <td>19.987330</td>\n",
       "      <td>-27.204305</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Flexible</th>\n",
       "      <td>-23.095745</td>\n",
       "      <td>19.954674</td>\n",
       "      <td>-27.323241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Adjusted (no cross-fit)</th>\n",
       "      <td>-27.137309</td>\n",
       "      <td>21.871141</td>\n",
       "      <td>-20.343291</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                     LATE       s.e.  % reduction\n",
       "Random Forest                  -19.591159  20.252981   -26.236779\n",
       "Gradient Boosting              -21.473448  21.536645   -21.561559\n",
       "Linear Regression              -27.272582  20.229335   -26.322902\n",
       "Lasso Baseline                 -22.005207  19.987330   -27.204305\n",
       "Lasso Flexible                 -23.095745  19.954674   -27.323241\n",
       "Linear Adjusted (no cross-fit) -27.137309  21.871141   -20.343291"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "med_est = np.median(np.array(estimates), axis=0)\n",
    "mean_se = np.array(std_err).mean(axis=0)\n",
    "adj_frame = pd.DataFrame(np.c_[med_est, mean_se], index=methods, columns=[\"LATE\", \"s.e.\"])\n",
    "adj_frame[\"% reduction\"] = (adj_frame[\"s.e.\"] - res_dataframe.loc[\"Food T_1\", \"s.e.\"]) * 100\n",
    "adj_frame[\"% reduction\"] /= res_dataframe.loc[\"Food T_1\", \"s.e.\"]\n",
    "adj_frame.loc[\"Linear Adjusted (no cross-fit)\"] = res_dataframe_adj.loc[\"Food T_1\"]\n",
    "adj_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0ioABDq3TK_o"
   },
   "source": [
    "## We now repeat the exercise for the other outcomes (excluding the simulation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "g-dLrCgrTm3n"
   },
   "source": [
    "Non-Food Consumption (Year 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "id": "zcKoNbJeTpSg"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>-3.262580</td>\n",
       "      <td>19.200077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient Boosting</th>\n",
       "      <td>-2.543165</td>\n",
       "      <td>20.288621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Regression</th>\n",
       "      <td>-5.154475</td>\n",
       "      <td>19.383861</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Baseline</th>\n",
       "      <td>-5.334120</td>\n",
       "      <td>19.471454</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Flexible</th>\n",
       "      <td>-6.341323</td>\n",
       "      <td>19.577609</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       LATE       s.e.\n",
       "Random Forest     -3.262580  19.200077\n",
       "Gradient Boosting -2.543165  20.288621\n",
       "Linear Regression -5.154475  19.383861\n",
       "Lasso Baseline    -5.334120  19.471454\n",
       "Lasso Flexible    -6.341323  19.577609"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Running Variable and Outcome\n",
    "investigated_outcome = \"conspcnonfood_t1\"\n",
    "df_ml = df.rename(columns={\"pov_index\": \"X\", investigated_outcome: \"Y\"})\n",
    "\n",
    "# Baseline covariates including consumption\n",
    "b_covs = df_ml.columns[[0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 18, 21]]\n",
    "\n",
    "# Fixed effects for localities\n",
    "i_fe = pd.get_dummies(df_ml['clus'], drop_first=True)\n",
    "\n",
    "# Flexible covariates including localities indicators\n",
    "f_covs = patsy.dmatrix('~ (' + ' + '.join(b_covs) + ')**2', data=df_ml, return_type='dataframe')\n",
    "\n",
    "# Dropping the intercept column that is automatically added by patsy\n",
    "f_covs = f_covs.iloc[:, 1:]\n",
    "\n",
    "Z_lasso = pd.concat([i_fe, f_covs], axis=1)\n",
    "\n",
    "h_fs = rdrobust(y=df_ml.Y, x=df_ml.X, masspoints=\"off\").bws.values[1, 0]\n",
    "\n",
    "\n",
    "def first_stage():\n",
    "    # Set up the cross-fitting\n",
    "    n = df_ml.shape[0]\n",
    "    Kf = 5  # Number of folds\n",
    "    sampleframe = np.repeat(range(0, Kf), math.ceil(n / Kf))\n",
    "    cfgroup = np.random.choice(sampleframe, size=n, replace=False)\n",
    "\n",
    "    # Matrix to store eta predictions\n",
    "    eta_fit = np.empty((n, 5))\n",
    "\n",
    "    # Create vector of observations to be considered in the first stage model\n",
    "    weights = np.abs(df_ml.X) < h_fs\n",
    "\n",
    "    for k in range(Kf):\n",
    "        fold = (cfgroup == k)\n",
    "\n",
    "        data_treated = df_ml.loc[(df_ml.X > 0) & ~fold & (weights > 0)]\n",
    "        data_control = df_ml.loc[(df_ml.X < 0) & ~fold & (weights > 0)]\n",
    "\n",
    "        data_fold = df_ml.loc[fold]\n",
    "\n",
    "        rf1 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        rf0 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 0] = (rf1.predict(data_fold[b_covs]) + rf0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lgbm1 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lgbm0 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 1] = (lgbm1.predict(data_fold[b_covs]) + lgbm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lm1 = LinearRegression()\n",
    "        lm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lm0 = LinearRegression()\n",
    "        lm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 2] = (lm1.predict(data_fold[b_covs]) + lm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        las_base1 = LassoCV()\n",
    "        las_base1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        las_base0 = LassoCV()\n",
    "        las_base0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 3] = (las_base1.predict(data_fold[b_covs]) + las_base0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        X_flex_treated = pd.concat([Z_lasso.loc[data_treated.index], data_treated[b_covs]], axis=1)\n",
    "        X_flex_control = pd.concat([Z_lasso.loc[data_control.index], data_control[b_covs]], axis=1)\n",
    "        X_flex_fold = pd.concat([Z_lasso.loc[data_fold.index], data_fold[b_covs]], axis=1)\n",
    "        X_flex_treated.columns = X_flex_treated.columns.astype(str)\n",
    "        X_flex_control.columns = X_flex_control.columns.astype(str)\n",
    "        X_flex_fold.columns = X_flex_fold.columns.astype(str)\n",
    "        las_flex1 = LassoCV()\n",
    "        las_flex1.fit(y=data_treated.Y, X=X_flex_treated)\n",
    "        las_flex0 = LassoCV()\n",
    "        las_flex0.fit(y=data_control.Y, X=X_flex_control)\n",
    "        eta_fit[fold, 4] = (las_flex1.predict(X_flex_fold) + las_flex0.predict(X_flex_fold)) / 2\n",
    "\n",
    "    return eta_fit\n",
    "\n",
    "\n",
    "eta_fit = first_stage()\n",
    "\n",
    "methods = [\"Random Forest\", \"Gradient Boosting\", \"Linear Regression\",\n",
    "           \"Lasso Baseline\", \"Lasso Flexible\"]\n",
    "\n",
    "\n",
    "def second_stage(eta_fit):\n",
    "    adj_results = []\n",
    "    for i in range(len(methods)):\n",
    "        M_Y = df_ml.Y - eta_fit[:, i]\n",
    "        rd_call = rdrobust(y=M_Y, x=df_ml.X, masspoints=\"off\")\n",
    "        adj_results.append([rd_call.coef.iloc[0].values[0],\n",
    "                            rd_call.se.iloc[2].values[0]])\n",
    "    return adj_results\n",
    "\n",
    "\n",
    "adj_frame = pd.DataFrame(second_stage(eta_fit), columns=[\"LATE\", \"s.e.\"],\n",
    "                         index=methods)\n",
    "adj_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4AR3BRv7Tp2D"
   },
   "source": [
    "Food Consumption (Year 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "id": "aoFEC4l-TsIV"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>62.169538</td>\n",
       "      <td>44.462430</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient Boosting</th>\n",
       "      <td>66.854620</td>\n",
       "      <td>45.057171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Regression</th>\n",
       "      <td>53.339545</td>\n",
       "      <td>45.131263</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Baseline</th>\n",
       "      <td>50.049116</td>\n",
       "      <td>45.546159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Flexible</th>\n",
       "      <td>54.287738</td>\n",
       "      <td>45.596947</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        LATE       s.e.\n",
       "Random Forest      62.169538  44.462430\n",
       "Gradient Boosting  66.854620  45.057171\n",
       "Linear Regression  53.339545  45.131263\n",
       "Lasso Baseline     50.049116  45.546159\n",
       "Lasso Flexible     54.287738  45.596947"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Running Variable and Outcome\n",
    "investigated_outcome = \"conspcfood_t2\"\n",
    "df_ml = df.rename(columns={\"pov_index\": \"X\", investigated_outcome: \"Y\"})\n",
    "\n",
    "# Baseline covariates including consumption\n",
    "b_covs = df_ml.columns[[0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 18, 21]]\n",
    "\n",
    "# Fixed effects for localities\n",
    "i_fe = pd.get_dummies(df_ml['clus'], drop_first=True)\n",
    "\n",
    "# Flexible covariates including localities indicators\n",
    "f_covs = patsy.dmatrix('~ (' + ' + '.join(b_covs) + ')**2', data=df_ml, return_type='dataframe')\n",
    "\n",
    "# Dropping the intercept column that is automatically added by patsy\n",
    "f_covs = f_covs.iloc[:, 1:]\n",
    "\n",
    "Z_lasso = pd.concat([i_fe, f_covs], axis=1)\n",
    "\n",
    "h_fs = rdrobust(y=df_ml.Y, x=df_ml.X, masspoints=\"off\").bws.values[1, 0]\n",
    "\n",
    "\n",
    "def first_stage():\n",
    "    # Set up the cross-fitting\n",
    "    n = df_ml.shape[0]\n",
    "    Kf = 5  # Number of folds\n",
    "    sampleframe = np.repeat(range(0, Kf), math.ceil(n / Kf))\n",
    "    cfgroup = np.random.choice(sampleframe, size=n, replace=False)\n",
    "\n",
    "    # Matrix to store eta predictions\n",
    "    eta_fit = np.empty((n, 5))\n",
    "\n",
    "    # Create vector of observations to be considered in the first stage model\n",
    "    weights = np.abs(df_ml.X) < h_fs\n",
    "\n",
    "    for k in range(Kf):\n",
    "        fold = (cfgroup == k)\n",
    "\n",
    "        data_treated = df_ml.loc[(df_ml.X > 0) & ~fold & (weights > 0)]\n",
    "        data_control = df_ml.loc[(df_ml.X < 0) & ~fold & (weights > 0)]\n",
    "\n",
    "        data_fold = df_ml.loc[fold]\n",
    "\n",
    "        rf1 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        rf0 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 0] = (rf1.predict(data_fold[b_covs]) + rf0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lgbm1 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lgbm0 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 1] = (lgbm1.predict(data_fold[b_covs]) + lgbm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lm1 = LinearRegression()\n",
    "        lm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lm0 = LinearRegression()\n",
    "        lm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 2] = (lm1.predict(data_fold[b_covs]) + lm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        las_base1 = LassoCV()\n",
    "        las_base1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        las_base0 = LassoCV()\n",
    "        las_base0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 3] = (las_base1.predict(data_fold[b_covs]) + las_base0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        X_flex_treated = pd.concat([Z_lasso.loc[data_treated.index], data_treated[b_covs]], axis=1)\n",
    "        X_flex_control = pd.concat([Z_lasso.loc[data_control.index], data_control[b_covs]], axis=1)\n",
    "        X_flex_fold = pd.concat([Z_lasso.loc[data_fold.index], data_fold[b_covs]], axis=1)\n",
    "        X_flex_treated.columns = X_flex_treated.columns.astype(str)\n",
    "        X_flex_control.columns = X_flex_control.columns.astype(str)\n",
    "        X_flex_fold.columns = X_flex_fold.columns.astype(str)\n",
    "        las_flex1 = LassoCV()\n",
    "        las_flex1.fit(y=data_treated.Y, X=X_flex_treated)\n",
    "        las_flex0 = LassoCV()\n",
    "        las_flex0.fit(y=data_control.Y, X=X_flex_control)\n",
    "        eta_fit[fold, 4] = (las_flex1.predict(X_flex_fold) + las_flex0.predict(X_flex_fold)) / 2\n",
    "\n",
    "    return eta_fit\n",
    "\n",
    "\n",
    "eta_fit = first_stage()\n",
    "\n",
    "methods = [\"Random Forest\", \"Gradient Boosting\", \"Linear Regression\",\n",
    "           \"Lasso Baseline\", \"Lasso Flexible\"]\n",
    "\n",
    "\n",
    "def second_stage(eta_fit):\n",
    "    adj_results = []\n",
    "    for i in range(len(methods)):\n",
    "        M_Y = df_ml.Y - eta_fit[:, i]\n",
    "        rd_call = rdrobust(y=M_Y, x=df_ml.X, masspoints=\"off\")\n",
    "        adj_results.append([rd_call.coef.iloc[0].values[0],\n",
    "                            rd_call.se.iloc[2].values[0]])\n",
    "    return adj_results\n",
    "\n",
    "\n",
    "adj_frame = pd.DataFrame(second_stage(eta_fit), columns=[\"LATE\", \"s.e.\"],\n",
    "                         index=methods)\n",
    "adj_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jGG6Nb_gTsdu"
   },
   "source": [
    "Non-Food Consumption (Year 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "id": "B9lCmaP1TvVz"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATE</th>\n",
       "      <th>s.e.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>47.496187</td>\n",
       "      <td>29.875396</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient Boosting</th>\n",
       "      <td>46.015727</td>\n",
       "      <td>30.694663</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Linear Regression</th>\n",
       "      <td>47.021684</td>\n",
       "      <td>29.760011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Baseline</th>\n",
       "      <td>46.433691</td>\n",
       "      <td>30.207670</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lasso Flexible</th>\n",
       "      <td>49.486267</td>\n",
       "      <td>30.457168</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        LATE       s.e.\n",
       "Random Forest      47.496187  29.875396\n",
       "Gradient Boosting  46.015727  30.694663\n",
       "Linear Regression  47.021684  29.760011\n",
       "Lasso Baseline     46.433691  30.207670\n",
       "Lasso Flexible     49.486267  30.457168"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Running Variable and Outcome\n",
    "investigated_outcome = \"conspcnonfood_t2\"\n",
    "df_ml = df.rename(columns={\"pov_index\": \"X\", investigated_outcome: \"Y\"})\n",
    "\n",
    "# Baseline covariates including consumption\n",
    "b_covs = df_ml.columns[[0, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 18, 21]]\n",
    "\n",
    "# Fixed effects for localities\n",
    "i_fe = pd.get_dummies(df_ml['clus'], drop_first=True)\n",
    "\n",
    "# Flexible covariates including localities indicators\n",
    "f_covs = patsy.dmatrix('~ (' + ' + '.join(b_covs) + ')**2', data=df_ml, return_type='dataframe')\n",
    "\n",
    "# Dropping the intercept column that is automatically added by patsy\n",
    "f_covs = f_covs.iloc[:, 1:]\n",
    "\n",
    "Z_lasso = pd.concat([i_fe, f_covs], axis=1)\n",
    "\n",
    "h_fs = rdrobust(y=df_ml.Y, x=df_ml.X, masspoints=\"off\").bws.values[1, 0]\n",
    "\n",
    "\n",
    "def first_stage():\n",
    "    # Set up the cross-fitting\n",
    "    n = df_ml.shape[0]\n",
    "    Kf = 5  # Number of folds\n",
    "    sampleframe = np.repeat(range(0, Kf), math.ceil(n / Kf))\n",
    "    cfgroup = np.random.choice(sampleframe, size=n, replace=False)\n",
    "\n",
    "    # Matrix to store eta predictions\n",
    "    eta_fit = np.empty((n, 5))\n",
    "\n",
    "    # Create vector of observations to be considered in the first stage model\n",
    "    weights = np.abs(df_ml.X) < h_fs\n",
    "\n",
    "    for k in range(Kf):\n",
    "        fold = (cfgroup == k)\n",
    "\n",
    "        data_treated = df_ml.loc[(df_ml.X > 0) & ~fold & (weights > 0)]\n",
    "        data_control = df_ml.loc[(df_ml.X < 0) & ~fold & (weights > 0)]\n",
    "\n",
    "        data_fold = df_ml.loc[fold]\n",
    "\n",
    "        rf1 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        rf0 = RandomForestRegressor(max_features=4, n_estimators=1000)\n",
    "        rf0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 0] = (rf1.predict(data_fold[b_covs]) + rf0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lgbm1 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lgbm0 = LGBMRegressor(verbosity=-1)\n",
    "        lgbm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 1] = (lgbm1.predict(data_fold[b_covs]) + lgbm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        lm1 = LinearRegression()\n",
    "        lm1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        lm0 = LinearRegression()\n",
    "        lm0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 2] = (lm1.predict(data_fold[b_covs]) + lm0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        las_base1 = LassoCV()\n",
    "        las_base1.fit(y=data_treated.Y, X=data_treated[b_covs])\n",
    "        las_base0 = LassoCV()\n",
    "        las_base0.fit(y=data_control.Y, X=data_control[b_covs])\n",
    "        eta_fit[fold, 3] = (las_base1.predict(data_fold[b_covs]) + las_base0.predict(data_fold[b_covs])) / 2\n",
    "\n",
    "        X_flex_treated = pd.concat([Z_lasso.loc[data_treated.index], data_treated[b_covs]], axis=1)\n",
    "        X_flex_control = pd.concat([Z_lasso.loc[data_control.index], data_control[b_covs]], axis=1)\n",
    "        X_flex_fold = pd.concat([Z_lasso.loc[data_fold.index], data_fold[b_covs]], axis=1)\n",
    "        X_flex_treated.columns = X_flex_treated.columns.astype(str)\n",
    "        X_flex_control.columns = X_flex_control.columns.astype(str)\n",
    "        X_flex_fold.columns = X_flex_fold.columns.astype(str)\n",
    "        las_flex1 = LassoCV()\n",
    "        las_flex1.fit(y=data_treated.Y, X=X_flex_treated)\n",
    "        las_flex0 = LassoCV()\n",
    "        las_flex0.fit(y=data_control.Y, X=X_flex_control)\n",
    "        eta_fit[fold, 4] = (las_flex1.predict(X_flex_fold) + las_flex0.predict(X_flex_fold)) / 2\n",
    "\n",
    "    return eta_fit\n",
    "\n",
    "\n",
    "eta_fit = first_stage()\n",
    "\n",
    "methods = [\"Random Forest\", \"Gradient Boosting\", \"Linear Regression\",\n",
    "           \"Lasso Baseline\", \"Lasso Flexible\"]\n",
    "\n",
    "\n",
    "def second_stage(eta_fit):\n",
    "    adj_results = []\n",
    "    for i in range(len(methods)):\n",
    "        M_Y = df_ml.Y - eta_fit[:, i]\n",
    "        rd_call = rdrobust(y=M_Y, x=df_ml.X, masspoints=\"off\")\n",
    "        adj_results.append([rd_call.coef.iloc[0].values[0],\n",
    "                            rd_call.se.iloc[2].values[0]])\n",
    "    return adj_results\n",
    "\n",
    "\n",
    "adj_frame = pd.DataFrame(second_stage(eta_fit), columns=[\"LATE\", \"s.e.\"],\n",
    "                         index=methods)\n",
    "adj_frame"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
